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ABSTRACT 

The thermal history of the intracluster medium (ICM) is complex. Heat input 
from cluster mergers, from AGN, and from winds in galaxies offsets and may even 
prevent the cooling of the ICM. Consequently, the processes that set the temperature 
and density structure of the ICM play a key role in determining how galaxies form. In 
this paper we focus on the heating of the intracluster medium during cluster mergers, 
with the eventual aim of incorporating this mechanism into semi-analytic models for 
galaxy formation. 

We generate and examine a suite of non-radiative hydrodynamic simulations of 
mergers in which the initial temperature and density structure of the systems are set 
using realistic scaling laws. Our collisions cover a range of mass ratios and impact 
parameters, and consider both systems composed entirely of gas (these reduce the 
physical processes involved), and systems comprising a realistic mixture of gas and 
dark matter. We find that the heating of the ICM can be understood relatively simply 
by considering evolution of the gas entropy during the mergers. The increase in this 
quantity in our simulations closely corresponds to that predicted from scaling relations 
based on the increase in cluster mass. 

We examine the physical processes that succeed in generating the entropy in order 
to understand why previous analytical approaches failed. We find that: (1) The energy 
that is thermalised during the collision greatly exceeds the kinetic energy available 
when the systems first touch. The smaller system penetrates deep into the gravitational 
potential before it is disrupted. (2) For systems with a large mass ratio, most of the 
energy is thermalised in the massive component. The heating of the smaller system is 
minor and its gas sinks to the centre of the final system. This contrasts with spherically 
symmetric analytical models in which accreted material is simply added to the outer 
radius of the system. (3) The bulk of the entropy generation occurs in two distinct 
episodes. The first episode occurs following the collision of the cores, when a large 
shock wave is generated that propagates outwards from the centre. This causes the 
combined system to expand rapidly and overshoot hydrostatic equilibrium. The second 
entropy generation episode occurs as this material is shock heated as it re-collapses. 
Both heating processes play an important role, contributing approximately equally 
to the final entropy. This revised model for entropy generation improves our physical 
understanding of cosmological gas simulations. 
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1 INTRODUCTION 

The structure of the X-ray emitting gas that is seen in galaxy 
groups and clusters is still not understood. This gas (the in- 
tracluster medium, hereafter ICM) dominates the baryonic 
mass content of clusters and is material which is left over 
from galaxy formation. As such its properties present an 
important clue to the galaxy formation process. If the pro- 
cesses that set its temperature and density structure could 
be understood they should provide valuable constraints on 
galaxy formation models. However, linking together models 
for galaxy formation and accurate numerical methods capa- 
ble of tracing the hydrodynamic evolution of the ICM has 
proved difficult. If cooling is neglected, the emission proper- 
ties of clusters can be calculated in a robust manner (e.g., 
Evrard et al. 1996); however such simulations cannot form 
galaxies and thus omit a fundamental component of the 
physics. On the other hand, if cooling is included in the sim- 
ulations, the results are not stable to changes in numerical 
resolution unless an effective form of heating is introduced 
to limit the cooling rate in small or early halos (e.g., Balogh 
et al. 2001; Borgani et al. 2006). Many different "feedback" 
mechanisms have been tried to oppose the cooling insta- 
bility. Examples include galaxy winds and thermal energy 
injection (e.g., Springel & Hernquist 2003), delayed cooling 
(e.g., Kauffmann et al. 1999), preheating of intergalactic gas 
(e.g., Kaiser 1991; Evrard & Henry 1991), thermal conduc- 
tion (e.g., Benson et al. 2003; Dolag et al. 2004), and heating 
by AGN (e.g., Dalla Vecchia et al. 2004; Sijacki & Springel 
2006). These schemes have all met with various levels of 
success, but have not yet achieved a simultaneous match to 
the observed galaxy luminosity function, stellar mass frac- 
tion and temperature and density structure of the ICM. One 
factor slowing the development of these theoretical models is 
the lack of understanding of the physical processes at work 
in the simulations. Because it is so difficult to quantify the 
likely effect of introducing new processes, the models must 
be developed on a largely trial and error basis. 

An alternative approach is to develop a semi-analytic 
model of the thermal history of the ICM, in which the 
complex physics is encapsulated in a small number of semi- 
empirical equations. This approach has been very successful 
in improving our understanding of galaxy formation (e.g., 
Cole et al 2000; Springel et al. 2005; Bower et al. 2006). 
The semi-analytic approach has been tried in several stud- 
ies (e.g., Wu et al. 2000; Bower et al. 2001; Benson et al. 
2003), but the results have been limited because the tech- 
niques for setting the gas distribution in the halos have been 
ad hoc. What is needed is to take a step back, and to better 
understand the results of simple gas hydrodynamic mod- 
els. A better physical understanding of how entropy of the 
ICM is set in these simulations would equip us with the 
techniques needed to attack the much greater physical com- 
plexity of the problem when cooling and galaxy formation 
are included. This is the subject of this paper. 

The starting point for the next generation of semi- 
analytical approaches to the ICM should be to describe 
the physical state of the gas in terms of its entropy dis- 
tribution function. Entropy is a powerful concept for un- 
derstanding the density and temperature structure of clus- 
ters, and was initially introduced as a means of quantifying 
the departures of cluster structures from the simple scal- 



ings expected for their mass and temperature distributions 
(Evrard & Henry 1991; Kaiser 1991; Bower 1997; Balogh et 
al. 1999; McCarthy et al. 2003). By knowing the entropy 
distribution of the gas, we can determine the density and 
temperature profiles of the ICM within a given dark matter 
gravitational potential (e.g., Voit et al. 2002). This process 
requires that an outer boundary condition is set (e.g. by 
computing the external pressure due to infalling gas) but the 
details of the boundary condition have only a weak effect on 
the profile. Discussing clusters in the language of entropy is 
extremely powerful, since the gas responds adiabatically to 
slow changes in the gravitational potential. Furthermore, the 
buoyancy of the ICM ensures that a relaxed system will have 
an ordered structure with the lowest entropy gas located in 
the deepest part of the gravitational potential. Only cooling 
and shock heating or mixing events alter the entropy distri- 
bution of the gas: cooling lowers the gas entropy, while shock 
heating and mixing can only raise it. Cooling in clusters is 
already well understood and can be modelled with simple 
semi-analytic techniques (e.g., McCarthy et al. 2004). How- 
ever, we need to understand how entropy is generated during 
cluster growth to self-consistently model the cosmological 
formation of structure in a semi-analytic way. 

So what sets the entropy distribution in clusters? Pre- 
vious papers have looked at how entropy is generated in 
smoothly infalling gas (Cavaliere et al. 1998; Abadi et al. 
2000; Tozzi & Norman 2001; Dos Santos & Dore 2002; Voit 
et al. 2003). For a spherically symmetric smooth shell, it is 
possible to determine the infall velocity and density of this 
material at the cluster virial radius. The bulk infall veloc- 
ity is converted into internal thermal energy at a shock sur- 
rounding the cluster. The entropy generated in the shock can 
be computed from the Rankine-Hugoniot shock jump condi- 
tions. Modelling the growth of clusters in this way produces 
realistic entropy profiles, however the normalisation of the 
entropy is somewhat too high, so that the model predicts 
clusters with average densities p gas , that are lower than ob- 
served (see Voit et al. 2003). 

However, the smooth accretion approximation is not 
valid in a cold dark matter (CDM) dominated universe 
since most of the mass accreted by clusters has previously 
collapsed into smaller virialised mass concentrations (e.g., 
Lacey & Cole 1993; Rowley et al. 2004; Cohn & White 2005). 
This leads to a problem for the picture above since the en- 
tropy generated drops rapidly if the accreted material is al- 
ready dense and warm prior to the accretion shock. Impor- 
tantly, if this deficit is present at every step in a system's 
merger history, the problem becomes greatly compounded 
over time. This effect is far larger than the overestimate of 
the cluster entropy obtained in the smooth accretion case. 
Modelling a realistic structure in the accreted material leads 
to systems that are much denser and more luminous in X- 
rays than observed systems. This issue is discussed by Voit 
et al. (2003). 

Recently, considerable progress has been made in sim- 
ulating the growth of clusters using purely numerical tech- 
niques (e.g., Borgani et al. 2004; Kravtsov et al. 2005). The 
numerical simulations clearly show that p gas « /;,ptot (where 
fb is the cosmic baryon fraction, and ptot is total baryon 
+ dark matter density) at large radii, in agreement with 
observations (e.g., Vikhlinin et al. 2006; McCarthy et al. 
2007). This result holds for a wide range of simulation pa- 



Shock Heating in Cluster Mergers 3 



rameters and is largely insensitive to numerical technique or 
resolution (Frenk et al. 1999). The agreement between sim- 
ulations and observations indicates that we do not yet have 
an adequate understanding of the lumpy accretion process, 
since analytic attempts to model this process yield results in 
strong discord with the observations. Our immediate aim is 
to improve our understanding of these simulations. Clearly, 
the approximate methods we wish to develop will not replace 
direct simulations, but they do provide an essential tool for 
understanding their results, for estimating the impact of lim- 
ited numerical resolution, and for exploring the vast param- 
eter space of heating mechanisms used in galaxy formation 
models. In subsequent papers, we will apply this understand- 
ing to improve the treatment of gas heating in semi-analytic 
codes galaxy formation models (e.g., Cole et al. 2000; Bower 
ct al. 2006). By introducing our models for the shock heat- 
ing of the ICM, we will be able to self-consistently model 
additional heating from supernovae and AGN. At present, 
the best efforts to track the impact of AGN heating in semi- 
analytic models (e.g., Bower et al. 2001) cannot follow prop- 
erly the development of entropy and use an ad hoc descrip- 
tion based on energy conservation. 

In this paper we address the problem of lumpy accre- 
tion by focusing on idealised simulations of merging sys- 
tems. We consider simple two-body collisions and explore 
how the entropy of the final system is generated. In setting 
up our two-body collisions, we remain as faithful as possi- 
ble to the results of cosmological simulations. We initialise 
our systems with structural properties guided by the results 
of recent cosmological simulations and explore a wide range 
of representative mass ratios and orbits. These simulations 
show that the spherical accretion model previously consid- 
ered is a poor representation of the true physical process in 
hierarchical cosmological models. We show that the kinetic 
energy of the infall lumps is much greater than previously 
estimated. Furthermore, we show that, if the mass ratio of 
the accretion event is large, only a small fraction of the ki- 
netic energy is dissipated in the smaller component. Most of 
the heating effect is felt by the larger system. Our simula- 
tions allow us to present a much improved physical model 
for entropy generation in clusters. 

The present paper is outlined as follows. In §2, we 
present the details of our simulation setup, including a dis- 
cussion of the initial structural conditions of our systems, 
the adopted orbital parameters of the mergers, and other 
characteristics of the simulations. In §3, we analyse the en- 
tropy history of the systems during the merging process and 
qualitatively compare it with the standard spherical accre- 
tion model. In §4, we present an analytic model that en- 
capsulates the essential physics of the merging process in 
our idealised simulations. Finally, in §5, we summarise and 
discuss our findings. 



2 SIMULATION SETUP AND METHODS 

The baryons in virialised systems formed in non-radiative 
cosmological simulations approximately trace the dark mat- 
ter with p gas ~ fbptot- Any reasonable model of shock heat- 
ing ought to reproduce this basic result. However, as out- 
lined in §1, it has so far proved quite difficult to construct a 
physical analytic model that can successfully pass this test. 



Exploration of idealised two-body mergers could provide the 
key insights necessary to achieve this goal. These idealised 
mergers should be representative of typical mergers in cos- 
mological simulations and should themselves reproduce the 
propagation of (near) self-similarity. However, this is by no 
means a guaranteed result. First, idealised approaches such 
as those adopted below may be too simplistic to mimic 
closely enough the typical merger event in cosmological sim- 
ulations. For example, the degree to which the self-similarity 
of systems formed in non-radiative cosmological simulations 
is sensitive to the exact initial conditions of the systems (i.e., 
structural properties) or to the properties of the merger or- 
bits is unclear. There are few studies in the literature to 
date that examine the dependence of the properties of the 
baryons on the detailed merger history of a system. On the 
other hand, the adopted resolution of our idealised mergers 
is typically better than that of cosmological simulations. So 
it is possible that we may resolve phenomena that are not 
yet adequately resolved in cosmological simulations. In ei- 
ther case, the reproduction of self-similarity in our idealised 
simulations is not an automatic result. Below, we explore 
under what conditions self-similarity is achieved in our ide- 
alised merger simulations. 

The following subsections present a description of our 
simulation setup and methods. This discussion goes into 
some detail and the reader who is mainly interested in the 
results of the simulations may wish to skip ahead to §3. 

2.1 Initial conditions 

We make use of the public version of the parallel TreeSPH 
code GADGET-2 (Springel 2005) for our merger simula- 
tions. The Lagrangian nature of this code makes it ideal 
for tracking the entropy evolution of specific sets of par- 
ticles (or in principle even on a particle by particle basis) 
over the course of the simulations. By default the code im- 
plements the entropy-conserving SPH scheme of Springel & 
Hernquist (2002). This is ideal for our purposes since the 
code explicitly guarantees that the entropy of a gas particle 
will be conserved during any adiabatic process. Thus, we 
are assured that any increase in this quantity is rooted in 
physics. 

We perform a series of binary merger simulations in- 
volving systems ranging in mass from 10 14 M Q < M200 < 
1O 15 M with the mass of the primary system (henceforth, 
we refer to the 'primary' system as the more massive of 
the two systems) in all collisions set by default to M200 = 
1O 15 M . (M200 is defined in §2.1.1.) Thus, we simulate colli- 
sions characterised by mass ratios ranging from 10:1 to 1:1. 
Although we have experimented with a variety of mass ra- 
tios in this range, we present the results for simulations with 
mass ratios of 10:1, 3:1, and 1:1 only, but note that these 
yield representative results. Even though we have chosen 
to focus on the high end of the mass spectrum, the results 
should also be applicable to mergers involving lower mass 
systems (e.g., galaxies) so long as the mass ratio is similar. 
This is by virtue of the fact that gravity, which is scale-free, 
is completely driving the evolution of the systems. One of 
the reasons for choosing to focus on clusters is that grav- 
itational processes dominate their final properties (X-ray 
observations demonstrate that non-gravitational processes 
such as cooling and AGN feedback influence the very cen- 
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tral regions only of massive clusters; see, e.g., McCarthy et 
al. 2007). So we expect that our idealised merger simula- 
tions will be directly applicable to the bulk of the baryons 
in real clusters, although we leave a comparison with X-ray 
(and/or SZ effect) observations for future work. As outlined 
in §1, our primary goal is to develop a physical description 
of entropy generation in mergers. 



2.1.1 Gas-only models 

We have set up and run collisions of systems composed 
purely of gas. Of course, observations demonstrate that dark 
matter dominates by mass the baryons in galaxies, groups, 
and clusters and is of fundamental importance in the process 
of structure formation. In this respect, it might be expected 
that the gas-only merger simulations will have limited ap- 
plicability. On the other hand, we anticipate that these ide- 
alised collisions will be considerably more straightforward to 
interpret than the gas + dark matter (hereafter, gas+DM) 
simulations (described below), given that the former have 
only a single phase to consider. If a physical understand- 
ing of the gas-only runs can be achieved then this could be 
quite helpful for understanding the role of dark matter in 
the combined runs. In fact, it is demonstrated in §3 that the 
gas+DM collisions exhibit what may be regarded as fairly 
minor deviations from the results of the gas-only runs. Thus, 
we find the gas-only simulations are potentially an excellent 
tool for elucidating the key physical processes in mergers of 
massive systems. 

The systems are constructed to initially be structural 
copies of one another. In particular, the gas is assumed to 
have a NFW density profile (Navarro et al. 1997): 



dP(r) 
dr 



GM (r) 



(3) 



p(r) = 



(r/r s )(l + r/r a ) 2 



where p a = M s /(47rrf) and 
M 200 



Ms 



ln(l + r 2 oo/r 3 ) - (r 2 oo/r s )/(l + r 2 oo/r a ) ' 



(1) 



(2) 



In the above, r 2 oo is the radius within which the mean 
density is 200 times the critical density, p cr it, and M200 = 
M(r 200 ) = (4/3)7^1,0 x 200p crit . 

To define the density profile for a halo of mass M200, 
we must set the scale radius, r„. The scale radius is often 
expressed in terms of a concentration parameter, C200 = 
r 2 oo/r a . There are numerous studies that have examined the 
trend between concentration parameter and mass in pure 
dark matter cosmological simulations (e.g., Eke et al. 2001). 
The fact that there is a trend at all implies that the dark 
matter halos in these simulations are not strictly self-similar. 
But the mass dependence of the concentration parameter is 
weak, varying by only a factor of 2 or so from individual 
galaxies to galaxy clusters (with lower mass halos being more 
concentrated). A fixed concentration parameter of C200 = 4, 
a value typical of galaxy clusters simulated in a ACDM con- 
cordance cosmology, is adopted for all of our systems. 

In order to completely specify the properties of the gas 
we must choose an entropy profile and an outer boundary 
condition. As outlined in §1, the gas entropy is probably 
the most useful quantity to track during the simulations. 
In order to specify the entropy profiles of the systems, it is 
assumed that the gas is initially in hydrostatic equilibrium: 



The entropjQ, K, is then deduced through the equa- 
tion of state P — Kp^. A boundary condition must be 
supplied before equation (3) can be solved and we specify a 
value for the pressure of the gas at r 2 oo, where we truncate 
the gas profiles (unless stated otherwise). There is some free- 
dom in our choice of the pressure at the edge of the system. 
However, for physically reasonable values of P(r2oo) there is 
in fact very little difference between the resulting profiles. 
For the bulk of the gas, we find that the requirement of 
hydrostatic equilibrium forces K(M sas ) into a near power- 
law distribution with an index of approximately 1.3. This 
can be understood by noting that at intermediate radii the 
NFW profile is approximated well by an isothermal profile 
(i.e., p oc r -2 ), which has an entropy profile characterised 
by d log K jd log M gas = 4/3. To ensure that our initial sys- 
tems are structural copies of one another and to ease the 
comparison between the initial and final systems later on, 
a value of P(r2oo) that establishes this near pure power-law 
entropy distribution all the way out to the system's edge is 
selected. In real systems, the hot diffuse baryons are sup- 
posedly confined by the ram pressure of infalling material. 
Experimenting with a physical model for this ram pressure, 
Voit et al. (2002) have demonstrated that groups and clus- 
ters are indeed expected to have near power-law entropy 
distributions out to large radii. 

The analytic profiles must be discretised into individual 
particles for input into the GADGET-2 code. In particular, 
the initial particle positions, velocities, and internal energies 
(per unit mass) must be specified. 

For the initial particle positions, we start by generat- 
ing a glass using the "MAKEGLASS" compiler option of 
GADGET-2. Basically, a Poisson distribution of particles is 
generated and is then run with GADGET-2 using — G in 
place of G for Newton's constant. We allow the simulation 
to run for ~1000 time steps to ensure the glass achieves a 
uniform distribution. This uniform distribution is then mor- 
phed into the mass profile that corresponds to the density 
profile given in equation (1). This is done by selecting a point 
inside the distribution of particles, which we have taken to 
be the centre of mass, ranking all of the particles according 
to the distance from this point, and then moving each par- 
ticle radially such that the desired mass profile is achieved. 

For the particle velocities, the baryons are initially as- 
sumed to be at rest in their hydrostatic configuration. 

The specific internal energy, /, defined as I — 
(3/2)P/p gas , of each particle is assigned by using the parti- 
cle's distance from the centre of the halo and interpolating 
within the analytic profiles derived above. 

Lastly, we surround the systems with a low density 
medium with a pressure set to P(r2oo). This medium is 
dynamically-negligible and is simply put in place to confine 
the gas particles near the system's edge. Without a confin- 
ing medium as much as 20% of the system's mass can leak 
beyond r 2 oo- 



1 We adopt this common re-definition of entropy, which is re- 
lated to the standard thermodynamic specific entropy via a log- 
arithm and an additive term that depends only on fundamental 
constants. 
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Figure 1. Testing the local Maxwellian approximation for the 
primary system. The dashed blue curve shows the initial mass 
profile extended out to r 2 5 . The solid red curve shows the resulting 
equilibrium mass profile after running the dark matter halo in 
isolation for 50 Gyr (i.e., many dynamical times). Over the range 
the 0.1 r/r200 1.0 the initial and final mass profiles follow 
each other closely. 



2.1.2 Gas + dark matter (gas+DM) models 

For the systems containing both gas and dark matter a 
slightly different procedure is used to construct the systems. 
First, we set up systems composed entirely of dark matter. 
The density profiles of the systems are given by equation 
(1). However, for reasons described below, we extend the 
dark matter well beyond r 2 oo, to an overdensity of 25 (for 
our adopted concentration r 2 5 ~ 2.44r 2 oo). As in the gas- 
only setup, a glass particle distribution is morphed into the 
desired mass profile. Unlike the gas, however, the dark mat- 
ter particles have no thermal pressure and must be assigned 
appropriate velocities in order to maintain this mass profile. 
To achieve this, we solve the the Jeans equation (see Binney 
& Tremaine 1987): 



d[ol{r)p{r)] _ GM(r) 



dr 



p(r) 



(4) 



for the velocity dispersion profile, oy(r), of the systems. 
Equation (4) implicitly assumes that the orbits of the dark 
matter particles are purely isotropic^. 

To solve equation (4) it is assumed that virial equilib- 
rium holds at the edge of the dark matter halo (see, e.g., 
Ricker & Sarazin 2001); i.e., 

2 Cosmological simulations demonstrate that pure isotropy is vi- 
olated in the outer regions of clusters. However, given that our 
(hydrostatic) gas-only simulations yield results that are remark- 
ably consistent with our gas+DM simulations (see §3.2), this im- 
plies that the resulting structural properties of our systems are 
not sensitive to how the initial (pre-merger) mass profiles are 
maintained. 



We have also experimented with the boundary condi- 
tion that Ur(r)p(r) — > as r — ► oo and obtained very simi- 
lar results within r 2 s. Henceforth, results are presented for 
the runs that make use of the boundary condition given by 
equation (5). 

Each of the three velocity components for a dark matter 
particle are assigned a random value picked from a Gaussian 
distribution with width equal to the 1-D velocity disper- 
sion, o>, at that radius. However, it is known that this local 
Maxwellian approximation does not result in a system that 
is initially in equilibrium (though it is not far removed from 
one). Of most concern is that a sharp truncation of the sys- 
tem at some finite radiu^f] results in a significant amount 
of the mass near this radius seeping out to much larger 
distance. Thus, the resulting equilibrium mass profile can 
be significantly different from that which was originally in- 
tended. In fact, this happens at all radii but in the system's 
interior the flux of particles moving to larger radii is offset 
by particles moving inwards (which were originally at larger 
radii). This obviously cannot happen near the system's edge 
where there are no particles originally at larger radii (and 
moving inwards) to replace those moving away from the 
system's centre. One way to overcome this problem is to 
apply a smooth tapering of the density profile beyond the 
system's edge and compute the distribution function of the 
particles as opposed to making the Maxwellian assumption 
(e.g., Kazantzidis et al. 2004; Poole et al. 2006). However, we 
adopt a more simplistic, but still effective, method for over- 
coming this problem. In particular, as mentioned above, the 
initial dark matter profile is extended well beyond r 2 oo, out 
to r 2 5 . The aim is to ensure that within r 2 oo , our region of in- 
terest, there is a large enough influx of particles to maintain 
the desired mass profile. 

The pure dark matter halos are run in isolation for many 
dynamical times. In Figure 1, we show that by extending 
the dark matter halo out to r 2 5, the resulting equilibrium 
mass profile for the primary system traces the initial (in- 
tended) profile within r 2 oo quite well. A slight deviation 
is apparent for r < 0.1r 2 oo, which is due to limited mass 
resolution. However, we explicitly demonstrate in the Ap- 
pendix that this deviation has negligible consequences for 
our equal mass mergers. In the case of our unequal mass 
mergers, the secondary systems are resolved with fewer par- 
ticles and therefore the deviation at small radii is slightly 
larger in these systems. However, as we demonstrate in §4.2 
and discuss further in the Appendix, the vast bulk of the en- 
ergy is thermalised in the more massive primary system. So, 
as long as the primary system is well resolved the properties 
of the final merged system should be robust. This is what 
likely accounts for the fact that the properties of massive 
virialised systems formed in non-radiative cosmological sim- 
ulations are robust to relatively large changes in numerical 
resolution (see, e.g., Frenk et al. 1999). 

A drawback of extending the halo out to much larger 
radii is that a significant fraction of the mass is located out- 
side the region of interest (thus, we potentially waste a good 

3 Note that some sort of truncation is necessary since the mass 
profile corresponding to equation (1) diverges at large radii. 
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deal of computational effort). In order to mitigate this issue, 
we take the equilibrium configuration after running the dark 
halos in isolation and simply clip particles that lie beyond 
rso. We have run the clipped halos for a further 10 Gyr and 
verify that the mass profiles within r2oo is virtually unaf- 
fected. This clipped equilibrium configuration is used to set 
up our initial gas+DM systems. We keep the current posi- 
tions and velocities for our initial values for the dark matter 
particles in the combined runs. For the gas particles we also 
use the positions of the dark matter particles (that is, for 
those particles that lie within r2oo) but reflect them through 
the centre of mass so that the gas particles do not lie directly 
on top of the dark matter particles. The gas particle veloc- 
ities are set to zero and the internal energy densities are 
determined by placing the gas in hydrostatic equilibrium 
within the total gravitational potential. As in the gas-only 
models described above, a pressure boundary condition is 
selected such that K(M ga _ s ) is nearly a pure power-law out 
to r 2 oo- 

The above procedure implies that we will have equal 
numbers of gas and dark matter particles in our combined 
runs within r2oo (i-e., double the number of particles used 
in the isolated dark matter run within that radius) . In order 
to conserve mass, so that the equilibrium state for the dark 
matter is still valid, we re-assign the dark matter particle 
masses to (1 — /&) times that used in the isolated run while 
the gas particles are assigned a mass of ft, times the dark 
matter particle mass in the isolated run. Here is the ratio 
of gas mass to total mass within r2oo, which is set to the 
universal value of Qt,/Sl m = 0.02/i~ 2 /0.3 = 0.136. 

The end result of the above procedure is that the gas 
traces the dark matter, both phases are initially in equilib- 
rium, and within r2oo both phases have a form that is very 
nearly the intended NFW distribution. As in the gas-only 
models, we surround the gas+DM systems with a low den- 
sity pressure-confining gaseous medium. 

2.2 Merger orbits and other simulation 
char acter ist ics 

We follow the approach of Poole et al. (2006) and use the 
results of cosmological simulations to help specify the or- 
bital properties of our idealised merger simulations. In par- 
ticular, we turn to the study of Benson (2005; but see also 
Tormen 1997; Vitvitskaet al. 2002; Wanget al. 2005; Khoch- 
far & Burkert 2006). Benson (2005) used a large collec- 
tion of N-body simulations carried out by the VIRGO Con- 
sortium to study the distribution of orbital parameters of 
substructure falling onto massive (primary) halos. Particu- 
lar emphasis was given to the 2-D distribution of the rela- 
tive radial and tangential velocities of the substructure as 
they passed through a spherical shell with a radius set to 
r — (1 ± 0.2)r v i r , where r v i r is the virial radius of the pri- 
mary halo. As anticipated, the peak of the distribution cor- 
responds to approximately the circular velocity of the pri- 
mary halo at the virial radius (see his Fig. 2). In particular, 
v = (v^ + f 2 ) 1 / 2 ~ l.lt; c (rvir), where v T and v t are the rela- 
tive radial and tangential velocities, respectively. For the in- 
dividual components, the peak of the distribution is approx- 
imately centred on v r ~ 0.9« c (?Vir) and vt ~ 0.65n c (r v i r ) 
and it is clear that there is a correlation between v r and Vt- 
However, as pointed out by Benson (2005), the centring of 



Table 1. Merger Orbital Parameters. 



Mass ratio 


Sim. type 


do 




ft 






(kpc) 


(km/s) 


(km/s) 


1:1 


eras-onlv 

£,£10 Will J 


4126 


1444 





3:1 


gas-only 


3494 


1444 





10:1 


gas-only 


3021 


1444 





1:1 


gas+DM 


4342 


1444 





1:1 


gas+DM 


4342 


1400.9 


350.2 


1:1 


gas+DM 


4342 


1291.6 


645.8 


3:1 


gas+DM 


3658 


1444 





3:1 


gas+DM 


3658 


1400.9 


350.2 


3:1 


gas+DM 


3658 


1291.6 


645.8 


10:1 


gas+DM 


3170 


1444 





10:1 


gas+DM 


3170 


1400.9 


350.2 


10:1 


gas+DM 


3170 


1291.6 


645.8 



the peak also depends on the total mass of systems (see his 
Fig. 4). In particular, the orbits of mergers involving mas- 
sive systems are more radial and less tangential than those 
of mergers involving low mass halos, presumably because 
more massive structures tend to form at the intersections of 
large-scale filaments. Unfortunately, the sample wasn't large 
enough to fully quantify this mass dependence. 

For the purposes of the present study, we adopt the fol- 
lowing set of orbital parameters, which are representative 
of the results of Benson (2005). Unless stated otherwise, 
a fixed relative total velocity corresponding to the circular 
velocity of the primary at T2oo, ~ 1444 km/s, is adopted ini- 
tially for all of our simulation^. This is slightly lower than 
the mean value found by Benson (2005), but consistent at 
the 1-sigma level. For the gas+DM simulations, we examine 
three different orbits for each mass ratio, corresponding to 
v t — 0, v t — v r /4:, and v t = Mr/2. We refer to these as the 
"head on" , "small impact parameter" and "large impact pa- 
rameter" runs, respectively. In practical terms, this implies 
v t ~ 0.243u c , p (r2oo) and v t ~ 0.447v c , p (r2oo) for the latter 
two. This choice spans the lower half of the vt/v r distribu- 
tion for massive halos seen in Fig. 4 of Benson (2005). For 
the gas-only simulations, which we run simply to help coax 
out the key physics during mergers, we simulate head on 
collisions only [i.e., v r = v c ,p{rwo) and vt = 0]. 

Each of the simulations are initialised with the gaseous 
components of the primary and secondary systems just 
barely touching. Thus, the initial separation, do, is just the 
summation of V2oo for the primary and r2oo for the sec- 
ondary. For the gas+DM simulations, this implies that there 
is initially some overlap of dark matter halos of the primary 
and secondary systems, which extends beyond the gaseous 
component. A complete summary of the adopted orbital pa- 
rameters of all of our simulations is presented in Table 1. 
Note that for a fixed mass ratio do differs slightly between 
the gas+DM and the gas-only simulations. This difference 
arises owing to the slight deviation of the mass profiles in 

4 We note that the circular velocity at r200 differs by only a small 
amount from the circular velocity at the virial radius, which is 
what Benson actually compares his velocities to. This just reflects 
the fact that the NFW mass profile gives rise to a nearly flat 
circular velocity profile at large radii. 
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the gas+DM simulations from the intended NFW distribu- 
tion. Finally, the simulations have been set up such that the 
physical and centre-of-mass (of the entire system) coordi- 
nates and velocities are identical, where we have approxi- 
mated the initial configuration as two point masses. 

We have carried out a mass resolution study of one of 
our simulations (see the Appendix). Based on this study 
we adopt the following conditions. The primary systems 
in our default runs have 50,000 gas particles within r2oo- 
This is the case for both the gas-only and gas+DM runs 
and implies the gas particle mass is 2 x 1O 1O M0 in the for- 
mer and f b (2 x 10 10 M Q ) = 0.272 x 10 10 M Q in the latter. 
In the gas+DM runs, the primary systems have ~76,500 
dark matter particles within rso and 50,000 within r2oo with 
a dark matter particle mass of (1 - f b )(2 x 10 10 Af Q ) = 
1.728 x 1O 1O M . These particle masses are held fixed for 
the secondary systems as well (thus, they contain fewer par- 
ticles than the primary). 

For the gravitational softening length, we adopt a fixed 
physical size of 10 kpc for both the gas and dark matter 
particles. (We have also experimented with softening lengths 
of 5 kpc and 20 kpc and found the differences in the entropy 
evolution to be negligible.) We use a fairly standard set of 
SPH parameters, including a viscosity parameter, avisc, of 

0. 8, a Courant coefficient of 0.1, and the number of SPH 
smoothing neighbours, iV S ph, is set to 50 ± 2. 

Finally, each simulation is run for a duration of 13 Gyr, 

1. e., approximately a Hubble time. Particle data is written 
out in regularly spaced intervals of 0.1 Gyr. 



3 SIMULATION RESULTS 

Below we present a detailed discussion of the entropy evolu- 
tion of our simulations. The reader who is mainly interested 
in a general understanding of this progression may wish to 
skip ahead to §3.3, which presents a summary of our find- 
ings. 



3.1 Gas-only simulations 

3.1.1 Entropy evolution 

We start by considering the evolution of the entropy in the 
gas-only merger simulations. Plotted in Figure 2 is the evo- 
lution for the 1:1 merger subdivided into spherical shells. 
Immediately apparent is the high degree of symmetry in 
the 1:1 merger. In all four regions the 25 th , 50 th , and 75 th 
percentiles follow each other quite closely. In addition, with 
the exception of the outermost ring, it is evident that the 
particles in both systems have achieved a nearly convergent 
state by the end of the simulation. Of course, it is possi- 
ble to continue running the simulation to determine exactly 
what the convergent state of the outermost regions will be. 
However, since we have already run the simulation for a 
Hubble time this implies that the outer regions of systems 
formed from similar mergers in cosmological simulations will 
also not have had sufficient time to completely virialise by 
the present day. Since one of our aims is to understand the 
results of such simulations, we limit the duration of our sim- 
ulations to 13 Gyr. Below we will compare this final config- 
uration with a scaled up copy of the initial systems. 




Figure 2. Entropy evolution in the 1:1 gas-only merger. The red 
and blue lines/hatched regions represent the primary and sec- 
ondary systems, respectively. The upper and lower bounds of the 
hatched regions represent the 75 th and 25 th percentiles, respec- 
tively, while the central (thick) lines represent the median. The 
panels show the evolution of the entropy of particles separated ac- 
cording to their initial position in the cumulative gas mass profile. 
Since the cumulative gas mass profile is a monotonic function of 
radius, the panels also separate the particles according to initial 
distance from the centre of their respective halos. 



Some discussion of how the particles in each system 
actually get their entropy is warranted. To help aid the dis- 
cussion we plot in Fig. 3 a sequence of snapshots of the 1:1 
merger with particles colour coded according to the frac- 
tional change in entropy since the previous simulation out- 
put. The general progression of the merger can be described 
as follows. Since the two systems are initially just barely 
touching and have a relative velocity of v CtP (r 200), they be- 
gin interacting as soon as the simulation starts. In particular, 
two large shock fronts are quickly established as the systems 
approach one another (see the panel corresponding to t = 1 
Gyr in Fig. 3). However, there is very little increase in the 
entropy of the particles in any of the four regions plotted 
in Fig. 2 until after the cores of the two systems collide at 
t ~ 2 Gyr into the simulation. The implication is that the 
initial shocks are actually quite weak (as the colour coding 
in Fig. 3 would also indicate). A plausible explanation for 
this behaviour comes when one compares the sound speed 
of the gas, c s , defined as 



(6) 



with the initial relative velocity of the merger. Since the ini- 
tial relative velocity of the merger is set to the circular ve- 
locity, we are effectively comparing the sound crossing time 
of a system with its dynamical time. The condition of hydro- 
static equilibrium ensures that the gas temperature will be 
such that these two timescales are comparable. In this case 
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Figure 3. Fractional change in particle entropy as a function of time for the 1:1 gas-only merger. Shown are particles in a slice through 
the centre of thickness 250 kpc (i.e., \z\ < 125 kpc). Particles are colour coded according to the fractional change in entropy since the 
last simulation output (0.1 Gyr ago). Black points are for a fractional change of less than 2%, blue are for a 2%-10% change, cyan are 
for a 10%-20% change, green are for a 20%-30% change, yellow are for a 30%-40% change, and red are for a > 40% change. For clarity 
the surrounding pressure-confining medium is not displayed. Each panel is 10 Mpc on a side. 



we find the mean sound speed of the systems is ~1200 km/s. 
Compared to the initial relative velocity (see Table 1), this 
implies an initial Mach number of only M w 1.2. Therefore, 
we should expect only mild shock heating to occur during 
the early stages of the merger (as observed in Fig. 2). How- 
ever, as we illustrate in §4, by the time the cores collide the 
relative velocity can approach or exceed nearly 3 times the 
initial value, effectively bringing the merger into the strong 
shock regime. 

Following the collision of the cores, a large shock wave 
is generated. This shock quickly propagates outwards, heat- 
ing the gas that was initially predominantly on the far side 
of each system and has not yet finished falling in (see the 
panels corresponding to t = 2 and 3 Gyr in Fig. 3). In fact, 
this shock wave, combined with the expansion of the shocked 
high pressure/density gas near the core of the merged sys- 



tem, succeeds in reversing the infall of the material. This 
outflow of gas proceeds nearly adiabatically for a period of 
time that is dictated by the initial distance of the particle 
from its system's centre. For example, for the outermost ring 
plotted in Fig. 2 this period of adiabatic expansion occurs 
from t « 2.5 Gyr for a duration of nearly 3.5 Gyr. Gas par- 
ticles initially located close to the centre of their respective 
system do not spend such a large time in this outflow phase, 
as they end up being more deeply embedded within in the 
merged system's potential well. 

The outflowing gas eventually halts and begins to re- 
accrete onto the core of the merged system. Although the gas 
accretes from virtually all directions, the symmetry of the 
1:1 merger is such that the accretion happens preferentially 
along and perpendicular to the original collisional axis. This 
occurs for a significant period of time, until the merged core 
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time (Gyr) time (Gyr) 



Figure 4. Entropy evolution in the 3:1 (left) and 10:1 (right) gas only mergers. The red and blue lines/hatched regions represent the 
primary and secondary systems, respectively. The upper and lower bounds of the hatched regions represent the 75 th and 25 th percentiles, 
respectively, while the central (thick) curves represent the median. The panels show the evolution of the entropy of particles separated 
according to their initial position in the cumulative gas mass profile. Since the cumulative gas mass profile is a monotonic function of 
radius, the panels also separate the particles according to initial distance from the centre of their respective halos. 



has achieved a near spherical symmetry. Following this, gas 
continues to accrete at the outskirts of the system but does 
so in a more or less spherical fashion. As evidenced from 
Figs. 2 and 3, the gas is slowly being shock heated during 
this period of re-accretion. The character of this episode of 
entropy production therefore differs significantly from the 
first abrupt episode following core collision. However, as Fig. 
2 indicates, both episodes are of comparable importance in 
terms of setting the final state of the gas. Finally, we note 
that the late time (nearly spherical) accretion of material 
is what gives rise physically to the continued increase of 
entropy in the outermost regions of the system until the end 
of the simulations (e.g., as in the bottom-right panel of Fig. 
2). 

An interesting question is whether or not much mixing 
occurs as a result of the merging process. We have tested 
this by tracking not only the entropy evolution of individ- 
ual particles but also their spatial distributions. Interest- 
ingly, even though the shock heating of the systems occurs 
in an "inside-out" fashion, we see very little evidence for 
mixing. For example, particles that were initially located 
near the centres of the two (pre-merger) halos (i.e., particles 
that initially have relatively low entropies) end up being lo- 
cated near the centre of the final merged system, whereas 
the large-radii (high entropy) particles end up occupying the 
outer regions of the final system. A plausible explanation for 
this behaviour is that the initial Mach number distribution 
of the particles varies relatively weakly with radius. This 
can be understood by considering the initial temperature 
profiles. The condition of hydrostatic equilibrium results in 
initial temperature profiles that vary by less than a factor of 
2 from the cluster centre to its periphery. This then implies 
the Mach number distribution varies by less than a factor of 



2 1//2 over the entire system. Therefore, to a large degree, one 
expects (and the simulations bear this out) that most of the 
particles will be shock heated to a similar degree. As a re- 
sult, convective mixing is minimal. This agrees well with the 
idealised cluster merger simulations of Poole et al. (2006) . 

The entropy evolution of the 3:1 and 10:1 gas-only merg- 
ers is plotted in Figures 4a and 4b, respectively. In a qual- 
itative sense, both the primary and secondary systems in 
these mergers behave in a similar manner to the systems in 
the 1:1 case. For example, they both generate two relatively 
weak shock fronts early on, with the secondary driving a 
bow shock into the primary and vice-versa. As in the 1:1 
case, they exhibit two main periods of entropy production, 
the first being associated with a strong shock approximately 
when their cores collide and the second being associated with 
an extended period of re-accretion shocks. Furthermore, we 
also see little evidence for mixing in the 3:1 and 10:1 mergers. 
Of course there are some differences between the three sets of 
simulations in detail. For example, because the sound speed 
of the gas in the secondary in the 10:1 collision is smaller 
than that in the 1:1 case, the initial Mach number for the 
secondary is higher. Consequently, the initial bout of shock 
heating before the cores collide is relatively more important 
for the secondary in high mass ratio mergers. 

It is clear from a comparison of Figs. 2 and 4a,b that 
the higher the mass ratio of the merger the more strongly 
heated the secondary is relative to the primary. This again 
may be tied to the fact that as one decreases the mass of the 
secondary its characteristic temperature decreases resulting 
in a decreased sound speed and an increased Mach num- 
ber. It is also expected if the preservation of self-similarity 
is achieved by heating both the secondary and the primary 
up to same level. To explain, K — P/p%s oc T/p ga s, but 
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Figure 5. Comparison of the resulting entropy distributions from the 1:1 gas-only merger with self-similar expectations. Left: Comparison 
of radial profiles. The short-dashed red curve represents the initial entropy profile of the primary, the long-dashed magenta curve represents 
a self-similarly scaled up version of this profile, and the solid green curve represents the final entropy profile of the total merged system. 
Right: Comparison of Lagrangian gas mass profiles. The short-dashed red, long-dashed blue, and solid green curves represent the final 
K(M gas ) distributions of the primary, secondary, and total merged systems, respectively. The horizontal dotted line represents the 
self-similar scaling. 



self-similarity means that all systems have the same inter- 
nal mass structure and therefore the characteristic entropy 
depends only on the temperature of the system. The virial 
theorem relates the temperature of a system with its mass 
via M oc T 3 / 2 , and therefore the characteristic entropy of a 
system scales simply as 



K oc M 



2/3 



(7) 



Initially, therefore, the ratio of the secondary's char- 
acteristic entropy to that of the primary is K s /K p — 
(M s /Mp) 2 / 3 . So the secondary requires extra heating to 
overcome its initial deficit compared to the primary if self- 
similarity is achieved by heating both systems to the same 
level. 



3.1.2 Final configurations 

How do the final entropy distributions of the gas-only merg- 
ers compare with self-similar expectations? Figure 5 presents 
this comparison for the 1:1 gas-only merger both in tra- 
ditional radial coordinates (left panel) and in physical gas 
mass coordinates (right panel). In Figure 5a, we compare 
the initial entropy radial profile of the primary system with 
the final radial profile of the total merged system. We also 
show a self-similarly scaled up version of the primary's initial 
entropy profile, achieved by scaling the entropy coordinate 
up by (Mtot/Mp) 2/3 = 2 2/3 and the radial coordinate by 
(Mtot/M p ) 1/3 = 2 1/3 . Fig. 5a shows that the scaled up ver- 
sion of the primary's initial entropy profile follows the final 
entropy profile of the merged system over approximately two 
decades in radius. This indicates that, by and large, the 1:1 
gas-only merger has scaled self-similarly through the merger. 



However, deviations are clearly visible, particularly at large 
radii. These are discussed in more detail below. 

The entropy distributions of observed systems are most 
commonly presented in radial coordinates (e.g., Piffaretti et 
al. 2005; Donahue et al. 2006; Pratt et al. 2006), and so 
it is useful to present theoretical models in such units if 
the goal is to explain the properties of observed systems. 
However, our immediate goal in this paper is to develop 
a physical shock heating model. For this purpose, the en- 
tropy distributions are best presented in Lagrangian (gas 
mass) coordinates. In Figure 5b, therefore, we plot the fi- 
nal -ft'(Mgas) distribution^] for the merged system and for 
the primary and secondary individually for the 1:1 gas-only 
merger. The gas mass coordinate has been normalised by 
the total gas mass in the systems (in this case 10 15 Mq for 
the primary and secondary and 2 x 1O 15 M0 for the final 
merged system). The entropy coordinate has been scaled to 
the self-similar expectations; i.e., the initial K(M gas ) distri- 
bution of the primary scaled up by a factor of (M to t /M p ) 2//3 . 
(Equivalently, we could use the initial entropy distribution 
of the secondary scaled up by a factor of [Mtot/Ms] 2 ^ 3 .) 
Perhaps the simplest way to achieve this state is by hav- 
ing the primary and secondary systems individually obey 
self-similarity following the merger. In this case, we would 
expect the final distribution of the particles belonging to the 
primary to be (Mtot/M p ) 2/ ' 3 times the initial distribution, 
while the secondary's final distribution would be a factor of 
(Mtot /Ms) 2//3 larger than its initial distribution. 



5 The physical meaning of the K(M gas ) distribution is most 
straightforwardly thought of in terms of its inverse, M gSlB (K), 
which is the total mass of gas with entropy lower than K. 
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Figure 6. Comparison of the resulting if(M gas ) distributions from the 3:1 (left panel) and 10:1 (right panel) gas-only merger with 
self-similar expectations. The short-dashed red, long-dashed blue, and solid green curves represent the primary, secondary, and total 
systems, respectively. The horizontal dotted line represents the self-similar scaling. 



If self-similarity is strictly obeyed, then the solid green 
curve in Figure 5b, which represents the final K(M gaa ) dis- 
tribution for the merged system, should lie on the horizontal 
dotted line. Fig. 5b shows that, without any fine tuning of 
the simulations, the central 80% of the gas mass of the final 
merged system lies within approximately 10% of the self- 
similar result. As expected, the symmetry of the 1:1 case is 
such that both the primary and secondary individually obey 
self-similarity as well. 

Beyond M gas /M gaSi tot ~ 0.8 or so there is an apparent 
entropy excess relative to the self-similar expectation. Be- 
low, in §3.1.3, we demonstrate that this effect is artificial 
and its origin can be attributed to the (unrealistic) abrupt 
truncation of the gaseous atmospheres of our idealised sys- 
tems at )"2oo- Fortunately, as demonstrated below, this effect 
is limited to the outermost regions of our systems only. 

The K(M gas ) distributions for the 3:1 and 10:1 gas-only 
mergers are plotted in Figures 6a and 6b, respectively. In 
both cases the final merged system is quite close to the self- 
similar result for the central ~ 70% of gas mass. Beyond this 
the edge effect discussed below kicks in. Thus, the 3:1 and 
10:1 cases both qualitatively and quantitatively resemble the 
1:1 case. Interestingly, however, the path through which the 
final merged systems in the 3:1 and 10:1 cases achieve self- 
similarity is by over-heating the primary and under-heating 
the secondary, not by both obeying self-similarity individ- 
ually. This implies that some of the infall energy initially 
associated with the secondary system went into thermalis- 
ing the gas of the primary system. We return to this point 
later. 

Finally, it is interesting in and of itself that the gas-only 
simulations obey self-similarity. This implies that the reason 
the baryons trace the dark matter in cosmological hydrody- 
namic simulations is not simply because the baryons are just 
following the orders of the gravitationally dominant dark 



matter. In the absence of dark matter the baryons would 
apparently still obey self-similarity. We revisit this result in 
§5.4. 

3.1.3 Excess Entropy at Large Radii 

The resulting K(M sas ) distributions of all our idealised 
merger simulations (including the gas+DM mergers dis- 
cussed below) show evidence for deviations from self- 
similarity beyond M gaa /M gaSit ot ~ 0.8 or so. Is this effect 
real or artificial? We hypothesise that this effect is artificial 
and is caused by the (unrealistic) abrupt truncation of our 
idealised systems at some finite radius. In particular, one ex- 
pects a SPH-based code to systematically underestimate the 
density of the gas particles near the edge of the truncated 
systems. This, in turn, will have the effect of overestimating 
the entropy boost these particles receive in shock heating 
events. 

To test the above, we try a simple modification to the 
1:1 gas-only merger simulation. In particular, we extend the 
initial gas profiles out to an overdensity of 100 (i.e., rioo) 
and re-run the simulation. The point is to see if initially ex- 
tending the gaseous atmospheres to larger radii will have the 
effect of quenching the (over-efficient) entropy production in 
particles near r2oo- To make a fair comparison between the 
extended and default runs, in setting up the extended gas 
systems we choose the pressure at rioo such that P(r2oo) is 
identical to that in our default run. This way we ensure that 
within r2oo the systems in both our default and extended 
runs are identical before being read into GADGET-2. An ad- 
ditional consideration is that of the infall velocity. Since the 
systems are more massive when extended (for our adopted 
concentration Mioo ~ I.26M200), the initial gravitational 
potential energy between the systems is larger (more nega- 
tive) than that for our default run (despite the fact that the 
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Figure 7. Testing the origin of excess entropy at large radii. 
Shown are the results of the 1:1 gas-only merger where the gaseous 
atmospheres have initially been truncated at r2oo (i- e -i the default 
setup) and at rioo- Extending the gaseous atmospheres to larger 
radii results in a more accurate density determination near T2oo 
which, in turn, mitigates the over-efficient entropy production 
seen at large radii in the default run. 



centres of the extended systems are initially separated by a 
larger distance; do — 2rioo). Thus, if the infall velocity when 
the centres are separated by 2r2oo is to be v c , P (r2oo) (as in 
the default run), a lower initial infall velocity is required. 
We calculate this velocity by assuming the systems may be 
approximated as point masses, which should be valid during 
the early stages of the merger. 

In Figure 7, we compare the final 7f(M ga s) distribu- 
tions for the merged systems (normalised to the self-similar 
result) in the default and extended runs. For the central 40% 
of gas mass, the K (M gas ) distributions are similar. However, 
beyond this point differences become readily apparent. In- 
terestingly, for the run where the initial gas profiles were ex- 
tended to rioo, the final A"(M gas ) distribution remains close 
to the self-similar result all the way out to M200. Beyond 
this point, a deviation from self-similarity is seen, as in the 
default run. However, since within M200 (or, equivalently, 
r2oo) the systems in the extended and default runs are ini- 
tially identical, we must conclude that the deviation from 
self-similarity seen at large in Fig. 5 (and similar figures) is 
indeed artificial and its origin is linked to a poor (SPH) es- 
timate of the true gas density near the edges of the idealised 
systems. 

Fortunately, the central regions of the default run are 
not significantly affected by how we treat the edge and there- 
fore it is this region that we focus on when comparing to 
self-similar expectations or our analytic models. In partic- 
ular, when constructing physical analytic models that de- 
scribe our simulations in §4, we restrict our focus to radii 
corresponding to M gas /M gaSit ot < 0.8. 
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Figure 8. Entropy evolution in the head on 1:1 gas+DM merger. 
The red and blue lines /hatched regions represent the primary and 
secondary systems, respectively. The upper and lower bounds of 
the hatched regions represent the 75 th and 25 th percentiles, re- 
spectively, while the central (thick) lines represent the median. 
The panels show the evolution of the entropy of particles sepa- 
rated according to their initial position in the cumulative gas mass 
profile. Since the cumulative gas mass profile is a monotonic func- 
tion of radius, the panels also separate the particles according to 
initial distance from the centre of their respective halos. 



3.2 Gas + dark matter (gas+DM) simulations 

3.2.1 Entropy evolution 

In Figure 8 we plot the evolution of the entropy with time 
for the head on 1:1 gas+DM simulation. Several differences 
are apparent between the evolution in this figure and that 
of the 1:1 gas-only simulation plotted in Figure 2. For ex- 
ample, the gas+DM run shows evidence for a small amount 
of entropy production right from the start of the simula- 
tion, whereas the gas-only run does not (at least not within 
the central regions). This difference may be attributed to 
the different methods used for constructing the initial con- 
ditions. In an azimuthally averaged sense, both the gas-only 
and gas+DM runs have, for our purposes, nearly identical 
initial conditions. However, because the initial gas particle 
positions in the gas+DM run were assigned based on the 
positions of dark matter particles (see §2.1.2), small local 
inhomogeneities are initially present in these runs. As a re- 
sult, the gas is not in perfect hydrostatic equilibrium and the 
systems undergo a slight re-adjustment when the simulation 
starts. Fortunately, the amount of entropy generated in this 
re-adjustment is small compared to the entropy produced in 
the merger shock heating and can be safely ignored. This is 
demonstrated below when we make an explicit comparison 
of the gas-only and gas+DM results. 

The most important changes between Figs. 2 and 8 re- 
late to the properties of the second (proper) period of en- 
tropy generation. In particular, the second phase of heating, 
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Figure 9. Fractional change in particle entropy as a function of time for the head on 1:1 gas+DM merger. Shown are particles in a slice 
through the centre of thickness 250 kpc (i.e., \z\ < 125 kpc). Particles are colour coded according to the fractional change in entropy 
since the last simulation output (0.1 Gyr ago). Black points are for a fractional change of less than 2%, blue are for a 2%-10% change, 
cyan are for a 10%-20% change, green are for a 20%-30% change, yellow are for a 30%-40% change, and red are for a > 40% change. For 
clarity the surrounding pressure-confining medium is not displayed. Each panel is 10 Mpc on a side. 



which was associated with a series of re-accretion shocks in 
the gas-only runs, begins earlier, rises more steeply, and con- 
tributes more to the final state of the system (particularly in 
the central regions) in the gas+DM run than in the gas-only 
run. What is the origin of this behaviour? 

To help answer this question, we plot in Fig. 9 a se- 
ries of snapshots colour coded according to the fractional 
change in entropy since the last simulation output (see Fig. 
3 for the analogous plot for the gas-only run). The evolu- 
tion is quite similar to that of the gas-only run early on. 
However, noticeable differences are present between the two 
at intermediate times (t ~ 4 — 5 Gyr). In particular, in the 
gas+DM run it is apparent that a second fast moving shock 
wave has been generated. This shock wave propagates out- 
ward in a nearly spherical fashion and resembles that of the 
first shock created when the cores collided. The origin of 



this second shock can be linked to the collisionless nature 
of the dark matter, which allows the dark matter cores to 
pass through one another, while the gas cores collide. Since 
the dark matter dominates the gas by mass, and the dark 
core is able to drag significant quantities of gas into the 
other hemisphere and away from the overall system's centre 
of mass. As a result, some of the gas undergoes a second 
period of collapse, collides with gas infalling from the other 
hemisphere, and produces the shocljf]. It is demonstrated 
later that the extra energy required for this second shock is 



6 In fact, the dark matter cores oscillate back and forth a number 
of times, each time dragging gas away from the overall system's 
centre and consequently generating more shocks. However, these 
additional shocks are extremely weak and generate virtually neg- 
ligible amounts of entropy. 
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derived from tapping the dark matter component (see, e.g., 
Fig. 21). Following this second shock, there is an extended 
period of re-accretion that proceeds in a similar fashion to 
that in the gas-only run. It is the combination of this second 
shock and the re-accretion phase that is responsible for the 
increased importance of the second major phase of entropy 
production in Fig. 8. 

The head on 3:1 and 10:1 gas+DM mergers also show 
the same qualitative trends. But the head on case is spe- 
cial and, therefore, the final state of such collisions may not 
be representative of typical mergers in cosmological simula- 
tions. On the other hand, cosmological simulations suggest 
that groups and clusters acquire the bulk of their mass by 
substructure falling in along filaments. So the head on sce- 
nario will not be as far removed from the typical merger 
event as in the case of collisions between galaxies. However, 
it is still important to quantify the differences between the 
head on and non-zero impact parameter cases. 

In Figures 10-15 we compare snapshots showing the spa- 
tial distributions of the primary and secondary particles and 
also the fractional change in entropy since the last simula- 
tion output for the all of the gas+DM merger simulations. 
We omit the panel showing the initial particle positions at 
t — 0, since they are the same for the three different orbital 
cases (thus, the energy to be thermalised is the same for all 
three cases). In these figures, the centre of each panel cor- 
responds to the overall centre of mass, with the secondary 
system initially approaching from the left and the primary 
from the right. In the non-zero impact parameter runs, the 
secondary's initial motion is towards the upper right, while 
the primary is initially moving towards the lower left. 

The general progressions of the 1:1 non-zero impact pa- 
rameter cases, presented in Figs. 10 and 11, are as follows. 
In the small impact parameter case (where v t = iv/4 ini- 
tially) the cores graze each other after t ~ 2 Gyr of infall. 
As a result, the cores are temporarily spared from signifi- 
cant shock heating. However, the systems are gravitation- 
ally bound and cannot avoid a major collision for long. In 
short order the cores cease moving apart and begin falling 
inward nearly radially, setting up an almost head on sec- 
ondary collision (at t ~ 4 Gyr) that greatly heats the core 
gas. Additional small shocks are generated by the back and 
forth sloshing of the dark matter cores, but in general their 
contribution to the entropy production is minor. As in the 
cases examined above, an extended period of re-accretion 
then ensues. In this case, the bulk of the accretion first hap- 
pens preferentially along and perpendicular to the axis of the 
secondary (near head on) collision. Qualitatively, the large 
impact parameter case proceeds much the same way. The 
main differences are as follows. The larger impact parame- 
ter means that the gas cores go virtually unheated during 
the first pericentric passage, as the cores of the two systems 
completely miss one another. The larger impact parame- 
ter also implies the amount of time spent between the first 
and second pericentric passages will be longer. Furthermore, 
the secondary collision is not directly head on in this case, 
meaning that some of the gas actually undergoes a third 
pericentric passage before settling. In both the small and 
large impact parameter cases, significant angular momen- 
tum is imparted to the gas and large, nearly circular, bulk 
motions remain present until the end of the simulations. 

Figs. 12-15 are completely analogous plots for the 3:1 



and 10:1 cases. In both head on collisions, the secondary 
essentially acts like a bullet, driving a large shock into the 
gas of the primary system and is able to easily penetrate 
all the way to the core of the primary. In fact, unlike the 
head on 1:1 case, the core of the secondary in the 3:1 and 
10:1 cases is able to remain somewhat intact even after a 
collision with the core of the primary. As a result, there 
is actually a period where the cores pass through one an- 
other. However, the tidal forces exerted on the secondary's 
core significantly stretch it along the collision axis. A second 
head on collision between the cores of the two systems then 
occurs. At the same time, the sloshing back and forth of 
the dark matter cores is generating small shocks. The small 
and large impact parameter 3:1 and 10:1 cases show evi- 
dence for even more complicated behaviour, with multiple 
collisions occurring and shock waves propagating in various 
directions simultaneously. 

3.2.2 Final configurations 

Given the wide variety of behaviours seen in Figs. 10-15, 
one might reasonably expect that the resulting gas proper- 
ties of the various mergers to be significantly different from 
one another. But an analysis of the final entropy distribu- 
tions indicates otherwise. Plotted in Figures 16a,b are the 
final A'(Mgas) distributions for the various 1:1 mergers. In 
Fig. 16a, it is demonstrated that the bulk of the gas in the 
primary, secondary, and total systems nearly maintain self- 
similarity. However, it is apparent that the small and large 
impact parameter cases show evidence for excess heating of 
the core gas relative to the head on case. This is illustrated 
most clearly in Fig. 16b, which compares the K(M sas ) distri- 
butions for the total systems of the three different gas+DM 
simulations. This plot shows that the cases with non-zero 
impact parameter heat the core at the expense of the outer 
regions of the systems. However, it is remarkable that only 
the central ~20% of the gas mass show significant devia- 
tions from the head on case given the large modifications to 
the initial orbital parameters. It is worth noting that non- 
radiative cosmological simulations also show evidence that 
the gas departs from a pure hydrostatic NFW profile in cores 
of systems (e.g., Frenk et al. 1999; Voit et al. 2005). So some 
deviation from self-similarity in our simulations might be ex- 
pected. 

Also plotted in Fig. 16b is the K (M gas ) distribution for 
the 1:1 gas-only merger. Outside the central regions, it traces 
the distribution of the head on gas+DM simulation remark- 
ably well. The two distributions begin to deviate from one 
another for M gas /Mg aSjt ot < 0.4 or so. The excess heating 
in the gas+DM can plausibly be attributed to energy ex- 
change between the gas and the dark matter (see, e.g., Lin 
et al. 2006). In §4, we show that typically 5-10% of the dark 
matter's energy is transferred to gas. This energy exchange 
occurs during the period when the dark matter cores of the 
primary and secondary are oscillating back and forth. 

In Figures 17a,b we compare the final K(M sas ) distri- 
butions for the various 3:1 mergers. Fig. 17a shows that the 
bulk of the gas for the total systems in all three orbit cases 
nearly obey self-similarity. As in the gas-only mergers dis- 
cussed in §3.1, self-similarity is achieved by over-heating 
the primary and under-heating the secondary. Similar to 
the 1:1 gas+DM mergers discussed immediately above, the 
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Figure 10. Fractional change in particle entropy as a function of time for the 1:1 gas+DM mergers. Top: Head on case. Middle: Small 
impact parameter case. Bottom: Large impact parameter case. Particle colour coding is the same as in Figure 9. Each panel is 6 Mpc on 
a side. 
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Figure 11. Spatial distributions of particles originally belonging to the primary (red) and secondary (blue) systems for the 1:1 gas+DM 
mergers. Top: Head on case. Middle: Small impact parameter case. Bottom: Large impact parameter case. For clarity, the secondary 
particles are plotted on top of the primary particles. Each panel is 6 Mpc on a side. 
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Figure 12. Fractional change in particle entropy as a function of time for the 3:1 gas+DM mergers. Top: Head on case. Middle: Small 
impact parameter case. Bottom: Large impact parameter case. Particle colour coding is the same as in Figure 9. Each panel is 6 Mpc on 
a side. 




Figure 13. Spatial distributions of particles originally belonging to the primary (red) and secondary (blue) systems for the 3:1 gas+DM 
mergers. Top: Head on case. Middle: Small impact parameter case. Bottom: Large impact parameter case. For clarity, the secondary 
particles are plotted on top of the primary particles. Each panel is 6 Mpc on a side. 
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Figure 14. Fractional change in particle entropy as a function of time for the 10:1 gas+DM mergers. Top: Head on case. Middle: Small 
impact parameter case. Bottom: Large impact parameter case. Particle colour coding is the same as in Figure 9. Each panel is 6 Mpc on 
a side. 
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Figure 15. Spatial distributions of particles originally belonging to the primary (red) and secondary (blue) systems for the 10:1 gas+DM 
mergers. Top: Head on case. Middle: Small impact parameter case. Bottom: Large impact parameter case. For clarity, the secondary 
particles are plotted on top of the primary particles. Each panel is 6 Mpc on a side. 
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Figure 16. The resulting K(M gSlB ) distributions for the 1:1 gas+DM mergers. Left: Comparison of the final entropy distributions for 
the primary (short-dashed red), secondary (long-dashed blue), and total (solid green) systems for the three different orbital cases. Right: 
Comparison of the final entropy distributions of the total systems in the gas+DM mergers for the head on (solid red), small impact 
parameter (short-dashed green), and large impact parameter (long-dashed blue) systems and the gas-only head on merger (dot-dashed 
magenta). The horizontal dotted line represents the self-similar scaling. 



3:1 mergers also show evidence for departures from self- 
similarity for the central 20% or so of the total gas mass. A 
comparison of the final distributions for the total systems in 
Fig. 17b illustrates this more clearly. Interestingly, even the 
head on 3:1 case shows evidence for strong central heating. 
In fact, the head on case appears to be even slightly more 
efficient at heating the core than the case characterised by 
a large impact parameter. (But we note that in the large 
impact parameter case there is still some residual heating 
occurring at the end of the simulation, so in the long run it 
may be the more efficient of the two.) A comparison of the 
entropy distributions of the head on gas+DM and gas-only 
mergers again highlights the fact that the two differ from 
each other only in the very central regions. 

Finally, in Figures 18a, b we show the same set of plots 
for the 10:1 gas+DM mergers. The trends in these plots 
follow those discussed above for the 1:1 and 3:1 mergers. 
The only difference that we would mention is that the 10:1 
cases exhibit a much lesser degree of central heating than 
the simulations discussed previously. 



3.3 Summary of Simulation Results 

In the next section we develop a physically-motivated ana- 
lytic model that attempts to encapsulate the key physics of 
the merging process just described. However, before present- 
ing this model we summarise the results of our simulations: 

• Both the primary and secondary systems gain the bulk 
of their final entropy through two distinct episodes. The 
first episode is associated with a strong, quickly propagating 
shock wave that is generated approximately when the cores 
of the two systems collide. The second episode of heating 



occurs over an extended period of time. In the gas-only sim- 
ulations, this second phase is driven entirely by a series of 
re-accretion shocks. In the gas+DM simulations, it is driven 
by a combination of re-accretion shocks and energy trans- 
fer from the dark matter to the gas, with the re-accretion 
shocks playing the dominant role. 

• In all of our simulations, the contributions of the first 
and second episodes of entropy production to the final dis- 
tribution are comparable. 

• We find that the bulk of the gas in our simulations 
matches the self-similar result to within 10%. Deviations 
from self-similarity are seen at large radii in both the gas- 
only and gas+DM simulations. This is almost certainly due 
to a truncation effect in our idealised setup (see §3.1.3). De- 
viations from self-similarity are also seen at small radii in 
the gas+DM simulations (which are real) and are at least 
partially due to energy exchange between the gas and dark 
matter. Some deviation might be expected at small radii, as 
the gas in systems formed in non-radiative cosmological sim- 
ulations also show departures there from a pure hydrostatic 
NFW distribution. 

• With the obvious exception of the symmetric 1:1 case, 
self-similarity of the final merged system is achieved by over- 
heating the gas in the primary while under-heating the gas 
in the secondary. This implies that some of the infall energy 
initially associated with the secondary has gone into heating 
the primary. 

Our simulations therefore differ markedly from the stan- 
dard spherical smooth accretion model where the ICM is 
built up one shell at a time from the inside out and where 
each shell is shocked a single time as it enters the virial 
radius. It seems remarkable that these two very physically 



Shock Heating in Cluster Mergers 19 




Figure 17. The resulting K(M gSlB ) distributions for the 3:1 gas+DM mergers. Left: Comparison of the final entropy distributions for 
the primary (short-dashed red), secondary (long-dashed blue), and total (solid green) systems for the three different orbital cases. Right: 
Comparison of the final entropy distributions of the total systems in the gas+DM mergers for the head on (solid red), small impact 
parameter (short-dashed green), and large impact parameter (long-dashed blue) systems and the gas-only head on merger (dot-dashed 
magenta). The horizontal dotted line represents the self-similar scaling. 
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Figure 18. The resulting K(M g!iB ) distributions for the 10:1 gas+DM mergers. Left: Comparison of the final entropy distributions for 
the primary (short-dashed red), secondary (long-dashed blue), and total (solid green) systems for the three different orbital cases. Right: 
Comparison of the final entropy distributions of the total systems in the gas+DM mergers for the head on (solid red), small impact 
parameter (short-dashed green), and large impact parameter (long-dashed blue) systems and the gas-only head on merger (dot-dashed 
magenta). The horizontal dotted line represents the self-similar scaling. 
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different scenarios yield structural properties that are as sim- 
ilar as they are. However, as discussed in detail by Voit et 
al. (2003), the spherical accretion model fails to match the 
results of cosmological simulations when it is modified to 
account for more realistic "lumpy" accretion. The reason 
for this failure is that increasing the density but keeping 
the total energy to be thermalised fixed results in decreased 
entropy production in the shock. However, our simulations 
show that significant shock heating does not occur until core 
collision. As a result, there is significantly more infall energy 
available for thermalisation. To a large extent, this boost 
offsets the otherwise reduced level of entropy owing to the 
increased density of infalling gas in the lumpy model relative 
to the smooth accretion model. 



DISCUSSION: 
RESULTS 



UNDERSTANDING THE 



4.1 A Simple Single-Shock Model 

We have demonstrated that our idealised merger simula- 
tions approximately preserve self-similarity, as seen in non- 
radiative cosmological simulations. As a result, we are now 
in a position to try to develop a physical analytic model for 
the entropy evolution of the primary and secondary systems 
in our idealised simulations. Taking our cue from the study 
of Voit et al. (2003), we consider a simple model whereby 
gas with some initial density p\ and entropy K\ is mov- 
ing with a velocity Vi n in the system's centre-of-mass frame. 
If this velocity is supersonic it will generate a shock front. 
Assuming that the post-shock gas is at rest in the centre- 
of-mass frame (i.e., that all of the energy associated with 
Vi n is thermalised), then the shock propagates with a ve- 
locity Wshock into the gas in the centre-of-mass frame. In 
the rest-frame of the shock, the gas therefore has a veloc- 
ity v\ — Vi n + v shock, and the post-shock gas has a veloc- 
ity V2 = "shock- The pre-shock (upstream) conditions are 
related to the post-shock (downstream) conditions via the 
well-known Rankine-Hugoniot jump conditions (e.g., Shu 
1992): 



p2 = vl _ (T+JjAji 

pi ~ v2 ~ (7 - 1)M? + 2 ' 

P 2 _ 2 7 (Mj - 1) 
Pi 7+1 
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(8) 
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where Mi =vi/c s is the Mach number. Equations (8) and 
(9) can be used to yield the jump condition relating pre- 
shock and post-shock entropy: 
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where we have used 7 = 5/3. Therefore, one can calculate 
the final entropy distribution if the Mach number of the 
shock is known. Unfortunately, it is non-trivial to measure 
the Mach number of a shock directly from the simulations. 
Shock heating is implemented in SPH simulations via an ar- 
tificial viscosity term which significantly broadens the shocks 
both spatially and temporally. This prevents one from easily 



applying equations (8)-(ll) to individual SPH particles (see 
Pfrommer et al. 2006). Another difficulty is that the Mach 
number is expressed in terms of pre-shock velocity in the 
rest-frame of the shock. So application of these equations 
to individual particles would require one to carefully track 
the evolution of the shock itself during the simulation. To 
avoid these difficulties, we use equation (8) to instead ex- 
press the Mach number in terms of the difference between 
the pre-shock and post-shock velocities: 
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Rearranging and replacing vi by Mic 3 , we obtain 

1 1/2 



Mi 



Mi 



Mi 



Mi 



5 „ 2/3 



(13) 



Therefore, given the initial entropy and density distri- 
butions of the gas and the velocity of the gas in the centre- 
of-mass frame, it is possible to solve the quadratic equation 
(13) for the Mach number of each particle. Deriving the 
final entropy distribution is then simply a matter of plug- 
ging these Mach numbers into equation (11). Of course in 
our idealised merger simulations we know precisely what the 
initial entropy and density profiles of the systems are, but 
what velocity do we pick for «i n ? Various possibilities exist, 
but Vin should not exceed the relative velocity of the cores as 
they are about to collide (i.e., when the gravitational energy 
between the two systems has been maximally converted to 
infall energy). In fact, the velocity will have to be quite a bit 
lower than this since, for example, in the 1:1 case using the 
maximum relative velocity for the gas in both systems would 
require twice as much energy as there is available to be ther- 
malised. Unfortunately, given the complicated nature of the 
simulations, the correct value of Vi n could fall any where be- 
tween zero and this upper bound. Previous analytic studies 
(e.g., Voit et al. 2003) adopted the infall velocity at the virial 
radius. However, it is evident from §3 that significant shock 
heating does not occur in our simulations until the cores of 
the two systems nearly collide. As a result, the infall velocity 
will be much larger than assumed in those analytic studies. 
Furthermore, it is also clear that energy is being exchanged 
between the primary and the secondary systems (for exam- 
ple, even in the 10:1 case the secondary is capable of driving 
a shock that significantly heats virtually all of the gas in 
the primary). This makes it even more difficult to assess a 
priori what is the appropriate value of Vi n to assign for the 
primary and secondary systems. 

We sidestep this problem by inverting the question: i.e., 
given the initial entropy and density profiles, what velocity 
is required to explain the final entropy distribution? To an- 
swer this question we simply try a range of different values 
for Vin and assess which gives the best match to the final 
entropy distribution. This velocity can be cast in terms of a 
requirement for the total amount of energy that must have 
been thermalised in order to explain the final entropy distri- 
bution. In §4.2, we examine whether or not there is enough 
bulk energy available to explain our findings. 

Although the simulated systems undergo two periods 
of entropy production, we start by trying to use the single- 
shock model outlined above to explain the observations. In 
principle this model, which just uses continuity and con- 
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Figure 19. Comparison of the final entropy distribution of the primary in the 1:1 (left) and 3:1 and 10:1 (right) gas-only simulation 
with our simple analytic shock heating model. The thick red curve represents the simulation result. The dashed curves represent the 
simple analytic model (equations 11 &i 13) for choices of v i n /v c , P (f"200 ) ranging from 0.75 — 2.00 in steps of 0.25 (bottom to top). 



servation equations to link upstream and downstream con- 
ditions, can effectively describe physical situations that are 
more complex than a single shock. Therefore it is the natural 
starting point for our investigation. 

We start by first examining the gas-only simulations, 
which make a useful benchmark for the more realistic 
gas+DM runs. In Figures 19a and 19b, the final entropy 
distributions of the primary and secondary systems in the 
1:1, 3:1, and 10:1 gas-only mergers are compared with the 
simple analytic shock heating model proposed above. For 
the analytic model, the upstream values of the entropy and 
density are taken from the initial conditions of the simula- 
tions (see §2). We try several different values of «m for both 
the primary and secondary systems, each corresponding to a 
unique prediction for the final entropy distribution of both 
systems. All curves plotted in Figs. 19a & 19b have been 
normalised to the self-similar expectation. In the discussion 
that follows, we exclude the high entropy tail at large values 
of M gas /M gaS: tot from consideration. The §3.1.3 shows that 
this tail results from the truncation of our idealised halos. 

Figs. 19a & 19b show that a centre-of-mass velocity 
ranging approximately from 0.9 < Vi n /v C:P (r2oo) < 1-25 is 
required to explain the final entropy distributions of the 
primary systems in the three runs. The secondary sys- 
tems require slightly lower velocities ranging from 0.75 < 
Vin/v c . p (r2oo) < 1.25. For both the primary and secondary 
systems there is a trend with the mass ratio of the merger, 
in the sense that the higher the mass ratio the lower the 
required velocity is to explain their final entropy distribu- 
tions. It is worth noting, however, that no single choice of 
Vin can explain the entire K(M sas ) profiles for the primary 
and secondary systems. In particular, if the analytic model 
is normalised to explain the intermediate regions of the en- 
tropy profiles (say 0.3 < M gas /M gaSi tot < 0.6), it system- 
atically under-predicts the level of the lowest entropy gas 



(Mgas/Mgas.tot < 0.2) compared to the simulations. Never- 
theless, the bulk of the gas in the primary and secondary 
systems can be adequately modelled by a fairly small range 
of velocities. 

An identical set of plots is presented in Figures 20a, b for 
the gas+DM simulations. Note that significantly higher ve- 
locities are required to explain the final entropy distributions 
of the primary and secondary systems in the gas+DM simu- 
lations. In particular, the primary systems typically require 
velocities ranging from 1.35 < Vi n /v CtP (r2oo) < 1.80, while 
the secondaries typically require 1.00 < Vin/v c ,p(»"20o) < 1-80 
(see Table 2) . Even though the systems in the gas+DM sim- 
ulations require significantly higher velocities than those in 
the gas-only simulations this does not necessarily imply that 
the energetic requirements of gas+DM mergers exceed those 
of the gas-only mergers. It should be kept in mind that for a 
system of mass M200 there is simply much more gas that re- 
quires heating in the gas-only simulations (the baryon to to- 
tal mass ratio of the gas-only simulations is ~7 times larger 
than the gas+DM simulations). Below we compare the en- 
ergy required by the simple shock heating model to match 
level of entropy production seen in the idealised simulations 
with our best estimates of the amount of energy available. 

4.2 Energy Considerations 

If we assume that the simple shock heating model provides 
a reasonable description of what is taking place in the sim- 
ulations, the velocity v ln can be used to estimate the total 
amount of energy that was thermalised in producing the fi- 
nal entropy distributions. Conveniently, v m is defined such 
that the post-shock gas is at rest in the centre-of-mass frame, 
so the total thermalised energy is just: 

E T ,tot = £ T ,p + E T ,s = "Mgas.p^n.p + ^^W^s (14) 
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Figure 20. Comparison of the final entropy distribution of the primary in the 1:1 (left) and 3:1 and 10:1 (right) gas+DM simulation 
with our simple analytic shock heating model. The thick red curve represents the simulation result. The dashed curves represent the 
simple analytic model (equations 11 & 13) for choices of v i n /v CtP (r2oo ) ranging from 0.75 — 2.00 in steps of 0.25 (bottom to top). 



Typically, this results in values ranging from approx- 
imately 2 — 6 x 10 64 ergs for the gas-only mergers and 

0. 5 — 2 x 10 64 ergs for the gas+DM mergers. So even 
though the gas+DM mergers require higher velocities to 
preserve self-similarity, their energy requirements are lower 
than those of the gas-only simulations. As mentioned above, 
this is simply because there is more gas that requires heating 
in the gas-only simulations. 

Interestingly, even though the energetic requirements 
are different between the two types of simulations, the way 
the energy is distributed does not appear to be. In partic- 
ular, in both types of simulations the thermalisation of the 
primary's gas dominates the thermalisation energy budget 
with 50%, «80%, and «95% of the total energy in the 1:1, 
3:1, and 10:1 mergers, respectively (see Table 2). The fol- 
lowing analytic relationship yields a remarkably good fit to 
our simulated mergers (including those involving mass ratios 
not presented in this paper): 

Interestingly, this is quite close to the case where the 
primary and secondary thermalise each other's infall energy, 

1. e., Et iP /Et, b ~ M p /M s . We also note that this type of 
scaling is naturally achieved if each gas particle thermalises 
the same fraction of the total available energy. This interest- 
ing result deserves further investigation, which we take up 
in the next paper of this series. 

An important consistency check of the simple shock 
heating model is to test whether or not there is enough en- 
ergy available to meet the requirements of the model. Hav- 
ing too much energy available isn't necessarily a problem, 
since there are numerous ways the available energy could be 
tapped (e.g., some could go into bulk kinetic circular mo- 



Table 2. Single-Shock Model Velocity/Energy Requirements 



Mp/M s 


Sim. type 


^in,p 


Vin,s 


Et,p/ Er,tot 






[fc, P (r 2 oo)] 


bc,p(?-200)] 




1:1 


gas-only 


1.25 


1.25 


0.50 


3:1 


gas-only 


1.25 


1.00 


0.82 


10:1 


gas-only 


0.9 


0.75 


0.94 


1:1 


gas+DM 


1.8 


1.8 


0.50 


3:1 


gas+DM 


1.65 


1.35 


0.82 


10:1 


gas+DM 


1.35 


1.00 


0.95 



tions of the gas or into the dark matter in the case of the 
gas+DM simulations). However, if there isn't enough energy 
available, the only possibility is that the model is incorrect 
or, at best, incomplete. 

With this in mind, we calculate the energy available to 
be thermalised in our idealised mergers. We do so via two 
different methods, with both yielding similar results. The 
first method, which we refer to as the "simulation method", 
takes advantage of the excellent energy conservation of the 
GADGET-2 simulations. The total energy of the gas is 

£tot,gas(t) = ^K gas (t) + £ Ugaa (t) + £i gaa (t) (16) 

where Ej^ sax (t), £"u gaa (i), and -Ei gaa (t) are the kinetic, po- 
tential, and internal (thermal) energies of the gas at time t. 
Stot, g as(t) is conserved in the gas-only simulations to better 
than 1%. Therefore, we can write 

£K gaa (t) + i?u gaB W + £i gafi W = i?K gaa o + £u gaa o + £i ga3 ,o(17) 

where £k s „,o, -Eu gaa ,o, and £i gaa ,o are the values of the three 
different energies at the start of the simulation. The total 
energy that is available to be thermalised at time t is just 
the initial (centre-of-mass) kinetic energy plus the change in 
the gravitational potential energy: 
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ET, tot (t) = E KgB . s ,o - [Ev g „(t) - -Eu gaa ,o] 



(18) 



which is relatively trivial to measure in the simulations. This 
energy estimate can be compared directly with the simple 
shock heating model estimate in equation (14). 

However, the above treatment is strictly only valid for 
the gas-only simulations. In the gas+DM simulations, en- 
ergy can be exchanged between the gas and the dark matter. 
But we can take advantage of the fact that the summation 
of the total energy of the gas and the total energy of the 
dark matter is conserved in these simulations. Thus, for the 
gas+DM simulations we modify equation (18) to read 

-E-r.tot = Ek s ^,o + [-Eti gas (t) - -Eu gas ,o] — E D M^ g ^(t) (19) 
where 



,(t) = [E KDM {t) + E VDM (t)] 
-[Ek dm ,o + Mj dm ,o] 



(20) 



is the energy exchanged between the gas and dark matter. 

In principle, this energy exchange can go either way, 
but in general we find that the dark matter loses energy to 
the gas. In Figure 21, we plot the evolution of the energy 
associated with the dark matter in the gas+DM simulations 
(see caption). By the end of the simulations, we find that 
the dark matter in the 1:1, 3:1, and 10:1 mergers has lost 
approximately 7-10%, 7-9%, and 2-3%, respectively, of its 
energy to the gas. Interestingly, this estimate doesn't de- 
pend much on the adopted orbital parameters, nor the mass 
resolution of the simulation. We will return to this issue of 
energy exchange below. 

The second method that we use to calculate the energy 
available to be thermalised, which we refer to as the "ana- 
lytic method" , is as follows. If we assume that the primary 
and secondary remain completely intact from their initial 
setup until the point when the centres of the two systems 
coincide, it is straightforward to calculate the total energy 
to be thermalised. Neglecting the thermal energy of the sys- 
tems (which remains fixed with time by construction), the 
total energy of the gas in the centre-of-mass frame is 



£tot,gas — c^^rel ^gas,ps(r"ps 



(21) 



where [i = M gaS;P ilf gas , B /(M gaS! p + M gas , s ) is the reduced 
mass, v re i is the initial relative velocity between the primary 
and secondary systems, r ps is the separation of their centres, 
and t/gas.ps is the gravitational potential energy. Assuming 
that the primary and secondary systems are rigid and ignor- 
ing the potential energy of the systems due to themselves 
(which, again, does not change by construction), the poten- 
tial energy of interaction between the two systems is 



£/gas,ps(r ps ) = — / </>pdMgas,s + - / <-.</.*/. ,.. ,, 



2 / Pgas.^pd'V + - / p gas , P <M V (22) 



where <f), the gravitational potential, is defined as 
p(x') 



(x) = -G 



(23) 
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Figure 21. Evolution of the total energy of the dark matter in 
the gas+DM simulations. The energy has been normalised to its 
initial value. Since the systems are bound the total energies are 
negative. Thus, a value of £ t ot,DM(*)/-E ; tot,DM(* = 0) greater 
than unity implies that energy has been lost to the gas, while 
a value less than unity means energy has been extracted from 
the gas. The green dot-dashed curve in the 1:1 panel is for the 
highest resolution simulation in the mass resolution study in the 
Appendix. 



Approximating the systems as point masses initially, 
when there is little or no overlap between them, the potential 
energy of the gas is just 



3 (r ps = d ) = -- — 



1 G_ 

2 do 
GfbMpM s 



+ Mp,gaaM s 



(24) 



do 



where the second line is true only if both the primary and 
secondary systems have the same baryon fraction. 
Therefore, the total energy of the gas is 

1 r , „, GfbMpM s 



E 



:p[Vc,p(r 2 oo)] 



do 



(25) 



The maximum energy available to be thermalised occurs 
when the centres of the primary and secondary coincide. 
We can therefore calculate the maximum energy available 
to be thermalised by subtracting the potential energy when 
the systems coincide (calculated by evaluating equation 22) 
from the total energy given in equation (25). 

In Figures 22a,b, we compare the simple shock heating 
model's energy requirements (see Table 2) with the energy 
available to be thermalised as estimated by both the simu- 
lation and analytic methods described above. We focus first 
on the gas-only results plotted in Fig. 22a. The solid curves 
represent the amount of energy available to be thermalised 
as estimated with the simulation method. The peak of the 
curves are reached at t m 1.5-1.8 Gyr, i.e., just slightly be- 
fore the cores of the primary and secondary collide. The 
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Figure 22. Comparison of the shock heating model's energy requirements with the maximum energy available to be thermalised in 
the gas-only (left) and gas+DM (right) simulations. In the panels on the right, the short-dashed red curves represent the results of the 
simulation method shifted up by Etot,DM(' = 13Gyr) — i?tot.DM(i = 0), i.e., by the total amount of energy lost by the dark mark to gas 
by the end of the simulation. For comparison, the dot-dashed blue lines represent the gas-only simulations, which have been renormalised 
by multiplying the energy by the baryon fraction of the gas+DM systems. 



amplitude of the peak is within a few tens of percent of the 
maximum energy estimated via the analytic method (dotted 
line) . This agreement indicates that our estimate of the max- 
imum energy available for thermalisation is robust and also 
demonstrates that one can estimate this energy reasonably 
well using simple analytic modelling. A comparison to the 
dashed line, which represents the energy requirement of the 
simple shock heating model described in §4.1, yields interest- 
ing results. Apparently, the 1:1 gas-only merger has sufficient 
energy available to accommodate the shock heating model's 
requirements. So the simple model provides a viable expla- 
nation for this collision. For the 3:1 gas-only merger there 
is a very small deficit of energy. However, in the case of the 
10:1 gas-only merger, both the simulation method and an- 
alytic method estimates of the maximum amount of energy 
there is to be thermalised fall short of the required amount 
by roughly a factor of two. 

Moving on to the gas+DM simulations in Fig. 22b, we 
find that there is insufficient energy available in any of the 
simulations to accommodate the requirements of the simple 
shock heating model (although the 1:1 is close). This is the 
case even when we account for the energy exchange between 
the gas and the dark matter. Thus, even though both the 
gas-only and gas+DM simulations preserve self-similarity, 
they do not give a consistent answer when compared to the 
simple analytic model. However, there are some similarities 
between the two in terms of their comparison with the an- 
alytic model. For example, both show a similar trend with 
mass ratio, in the sense that agreement gets worse for higher 
mass ratios. Furthermore, the energy shortfall is less than 
about a factor of 3 for all of the simulations we have per- 
formed. Thus, while the simple analytic model fails to ex- 



plain the results, it does not fail by a huge margin. This 
motivates us to consider modifications of the model. 



4.3 A Double-Shock Model 

In §3 we found that there are actually two major periods 
of entropy production experienced by the primary and sec- 
ondary systems in our merger simulations. We now seek to 
modify the simple shock heating model presented in §4.1 to 
account for this behaviour. The relevant question is, for a 
fixed amount of energy to be thermalised, is a double-shock 
model capable of generating more entropy than a single- 
shock model? In the case where the post-shock conditions, 
as dictated by the jump conditions, simply become the pre- 
shock conditions for the second shock, we find that the an- 
swer is 'no'. In general, we find that as one increases the 
number of shocks over which the energy is to be thermalised, 
the resulting final entropy decreases. 

However, it quickly becomes apparent from an examina- 
tion of the simulations that the properties of the gas evolve 
significantly between the end of the first shock and the onset 
of the second. In particular, there is a period of adiabatic 
expansion between the two shocks which likely arises as a 
result of the fact that not all the kinetic energy was ther- 
malised in the first shoclfl The net result is that the typical 
density of the gas is significantly reduced between the shocks 

7 The expansion could also be partially due to a re-adjustment 
of the gas towards a new hydrostatic configuration. Note, how- 
ever, that this cannot be the whole story since no further shock 
heating would be expected in this case. A period of re-accretion 
is required. 
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Figure 23. Evolution of the median density for particles that 
were located within the spherical shell 0.25 $C M gas /M gaSi tot Ss 
0.75 initially. The solid red, dashed green, and dot-dashed blue 
curves represent the 1:1, 3:1, and 10:1 mergers, respectively. For 
the gas+DM simulations only the head on results are plotted. 
Note for the primary systems that the density drops below its 
initial (pre-merger) value between the end of the first shock and 
the onset of the second shock at t ~ 4 Gyr. For the secondary sys- 
tems the density increases (except for the 1:1 case), as expected. 
However, only a small fraction of the total energy is thermalised 
in the secondary. 

and can even drop below its pre-merger value (see Fig. 23). 
Furthermore, the drop is largest for the highest mass ratio 
mergers, precisely where we find the largest energy deficits 
between the simulations and the single-shock model. Drop- 
ping the density between the shocks will have the effect of 
increasing the amount of entropy generated in the second 
shock relative to the case where there is no adiabatic expan- 
sion between the shocks. More importantly, is the decrease 
in density between shocks large enough to generate more 
entropy than the single-shock model? 

To answer this question, we have compared the single- 
and double-shock models head to head for an idealised parcel 
of gas with an initial density pi n j t and an initial entropy 
^init- In particular, in Fig. 24 we plot three sets of curves 
representing three different comparisons, each characterised 
by a different total amount of energy to be thermalised. The 
total energies have been chosen such that, in the single-shock 
model, the (shock frame) Mach numbers are 1.5, 2.5, and 
3.5 (bottom to top). For the single-shock model, the Mach 
number is all that is required to predict the ratio of final to 
initial entropy (see equation 11). The three horizontal dotted 
lines in Fig. 24 represent the predictions of the single-shock 
model. 

For the double-shock model, we use the entropy and 
density evolution plots of the simulations as guides. Exam- 
ination of the evolution of the entropy of our mergers (e.g., 
see Figs. 2, 4, and 8) indicates that the first and second 
shocks contribute comparable amounts of entropy to the fi- 
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Figure 24. A comparison of the entropy generated by the single- 
and double-shock heating models. Horizontal dotted lines show 
the predictions of the single-shock heating model for (shock 
frame) Mach numbers of 1.5, 2.5., and 3.5 (bottom to top). The 
solid curves show the predictions of the double-shock model as 
a function of the density between the two shocks. Note that if 
the density drops below fs 40% of its pre-merger value then the 
double-shock model generates more entropy than the single-shock 
model for a fixed amount of energy to be thermalised. 



nal state. To first order, therefore, we surmise that the first 
and second shocks thermalise comparable amounts of energy 
(i.e., for the purposes of this toy model, we assume each 
shock thermalises half of the total energy.) Between the first 
and second shocks, there is a period of adiabatic expansion 
during which the density drops to some value p\^,i. Using 
the new density and entropy, we can solve for thermal en- 
ergy of the gas (or, equivalently the sound speed). The ratio 
of the thermal energy to the remaining energy to be ther- 
malised (i.e., half the initial energy) sets the Mach number 
of the second shock which, in turn, allows us to compute the 
final entropy. 

In Fig. 24, the solid curves show the predicted trend 
between final entropy and density between the two shocks for 
double-shock model. A comparison between the single- and 
double-shock models demonstrates that if the density drops 
below about 40% of its initial pre-merger value then the 
double-shock model does indeed generate more entropy than 
the single-shock model. This is quite promising since Fig. 
23 demonstrates that the primary systems, which dominate 
the thermalisation budget, have their densities reduced to 
at least this level (and lower for the 3:1 and 10:1 cases). 

Given these results, we test the model further by tai- 
loring the total energy to be thermalised (and, therefore, 
the Mach numbers) to match our merger simulations more 
closely. In particular, we assume the total energy is taken to 
be the peak of the solid curves plotted in Figs. 22a,b corre- 
sponding to the "simulation method" estimate. For simplic- 
ity, we further assume that each gas particle thermalises the 
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Figure 25. Final entropy predicted by the double-shock model as a function of the density pi^2 between the two shocks for the gas-only 
(left) and gas+DM (right) simulations. The solid red, dashed green, and dot-dashed blue curves represent the 1:1, 3:1, and 10:1 mergers, 
respectively. The entropy has been normalised to the expected self-similar result, while the density has been normalised to its initial 
pre-merger value. The vertical dotted lines show the minimum density estimated from Fig. 23 at t m 4 Gyr. 



same amount of energy. So, for example, the primary sys- 
tem thermalises M P /[M P + M s ] times the total energy. As 
above, the Mach numbers are set by computing the ratios 
of thermal energy to energy to be thermalised for the first 
and second shocks. 

In Figures 25a,b we plot the final entropy predicted by 
the double-shock model for the primary systems in our simu- 
lations. The entropy has been scaled to the self-similar result 
while pi->2, the density between the two shocks, has been 
scaled to the initial pre-merger density. We use the vertical 
dotted lines, which are meant to represent the density min- 
imum between the shocks seen in Fig. 23 (at t » 4 Gyr), to 
select the appropriate predicted entropy. 

A comparison of the predicted entropy in Figs. 25a,b 
with that of the primary in our simulations in Figs. 5-6 (gas- 
only) and Figs. 16a, 17a, and 18a (gas+DM) demonstrates 
that in all cases, including the 10:1 mergers, the double- 
shock model can reproduce the simulation results to within 
10%. Furthermore, we have found that this result is not very 
sensitive to the way in which we've distributed the energy 
over the two shocks. For example, we obtain very similar 
results if the second shock thermalises anywhere between 
20-80% of the total energy, as opposed to half. 

The double-shock heating model therefore appears to 
provide a simple framework for understanding the typical 
level of entropy generated in our mergers. In §5.3, we present 
a point by point algorithm that the reader can use to apply 
the double-shock model. 



5 CONCLUSIONS 

We have presented a series of simulations aimed at exploring 
the generation of entropy during cluster mergers. The results 



show that the entropy generated is remarkably robust. We 
show that the generation of entropy is largely independent 
of the impact parameter of the collision, and that similar 
results are obtained for simulations only involving gas and 
for mergers of systems that contain a cosmological mixture 
of gas and dark matter. The resulting entropy profiles also 
depend little on the resolution of the simulation, once more 
than ~ 10 4 particles are placed in each halo. These results 
hint that a general principle is at work, and that the genera- 
tion of entropy can be understood as a general process that 
converts the gravitational potential energy released during 
the collapse of the system into the thermal energy of the 
ICM. Based on this reasoning, it should be possible to de- 
velop a simple model for the evolution of the entropy distri- 
bution as clusters grow in mass through mergers by exam- 
ining the potential energy released in the collapse and the 
efficiency with which this is converted into thermal energy. 

5.1 Equal mass mergers 

We explore mergers of systems with a variety of mass ra- 
tios. In each case, we find that the entropy generated is ap- 
proximately sufficient to place the final system on the same 
entropy scaling relation that was used to generate the orig- 
inal system. During equal mass mergers symmetry ensures 
that both systems are heated to the same degree. For the 
bulk of the gas, the powerlaw slope of the resulting entropy 
profile changes little compared to the original (see, e.g., Fig. 
5a) and the degree of heating raises the normalisation of the 
entropy profile by » 2 2 / 3 . As a result the final entropy dis- 
tribution is scaled so that the final system is a self-similar 
copy of the original. 

The exact process by which this entropy is generated is 
far from simple. The two clusters are in contact with each 
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other at the start of the simulation, however, their infall ve- 
locity is not sufficiently high compared to the sound speed 
to generate a strong shock. The system becomes highly com- 
pact along the infalling axis with compressed material tend- 
ing to flow out along the orthogonal plane. A strong shock is 
not generated until the cores of the two systems are super- 
posed. At this point, a strong shock is generated, propagat- 
ing rapidly out from the central regions. However, not all of 
the available infall energy is thermalised in this first shock. 
Some of it remains in kinetic form and succeeds in driving 
a period of adiabatic expansion. Eventually, however, the 
remaining energy is thermalised in a series of shocks as ma- 
terial is re-accreted by the merger remnant. 

The generation of entropy thus shows two distinct 
peaks. The first corresponds to the strong shock generated 
as the cores collide, the second corresponding to the re- 
accretion of material that tried initially to escape from the 
system. 

The delay to the initial shock in this system plays an 
important role in determining the entropy generated. Be- 
cause the system collapses prior to the first shock, consider- 
ably greater binding energy is available to be thermalised. 
Superposing the initial mass distributions provides a good 
estimate of the available binding energy, and modelling this 
energy as being thermalised in a single strong shock provides 
a reasonable approximation to the entropy generated in the 
equal mass mergers. 

5.2 Unequal mass mergers 

For mergers between unequal mass halos, the generation of 
entropy is distributed unequally. Visually, we see that the 
smaller component remains essentially intact as it plunges 
into the centre of the main system (see, e.g., Figs. f3 and 15). 
We find that, although the kinetic energy of the collapsing 
system is primarily localised in the smaller mass system, this 
energy is largely thermalised in the more massive progeni- 
tor. As a result, the heating of the more massive progenitor 
exceeds what is predicted by the self-similar scaling relation 

K oc M 2/3 , 

while the heating of the less massive component falls short 
of that needed for self-similarity. Despite this, the over- and 
under-heating of the components combine is such a way that 
the final system comes close to following the self-similar re- 
lation. We provide an analytic fit to the ratio of the energy 
dissipated in the two components in equation (15). This is 
close to assuming that the energy is exchanged between the 
primary and secondary components. 

Thus we find that the mergers tend to produce scaled 
up copies of the original systems. This is good news since 
we started with systems having properties close to those of 
observed clusters, and we chose cosmologically likely values 
for the infall velocities. Our simulations show that the nor- 
malisation of this entropy profile is not a coincidence. Given 
the infall velocity distribution expected in a CDM universe 
this profile is a stable configuration. 

As with the 1:1 merger case, we can estimate the energy 
that is available to be thermalised by tracking the evolution 
of the potential energy. This energy significantly exceeds the 
initial kinetic energy of initial system, showing that the long 



survival time of the secondary is responsible for much of the 
entropy generation. However, we find that a single shock 
model for the thermalisation of this energy under-predicts 
the entropy generated, particularly in the case of the 10:1 
mass ratio. In order to match the entropy generation we find 
that it is necessary to model the two shock process that is 
seen in the simulations. A key ingredient of the entropy gen- 
eration is the drop in the gas density as the system responds 
to the first shock. As the remaining binding energy is ther- 
malised in this more diffuse medium, the entropy generation 
is more efficient. 

5.3 An algorithm for computing shock heating 

A good way to summarise the findings of this paper is 
to sketch out an algorithm for computing the entropy 
generated during the merger. Future papers will present the 
results from implementing this within cosmological merger 
trees. We summarise the algorithm as follows. 

1. ) Calculate the energy available for thermalisation. 
Cosmological simulations suggest that the secondary (less 
massive) system crosses the virial radius of the primary 
system with a total (relative) velocity of approximately the 
circular velocity of the primary system at that radius. The 
total energy is therefore given by equation (25). (Note, a 
similar energy estimate can be derived by calculating the 
potential energy between the systems at turnaround, when 
the relative kinetic energy is zero.) The energy available for 
thermalisation can then be obtained by subtracting from 
this the potential energy between the systems when their 
centres coincide (see §4.2). 

2. ) Distribute this energy in appropriate proportions 
to the primary and secondary systems. Our simulations 
indicate that the bulk of the energy is thermalised in the 
more massive primary system. Equation (15) provides a fit 
to how the energy should be divided up as a function of 
mass ratio. 

3. ) If the merger is 1:1 (or very nearly so), calculation of the 
post-shock properties is now straightforward. The energy 
estimated above can be converted into an estimate for 
the (centre-of-mass) velocity, v in , for both systems. This, 
in turn, may be converted into an estimate of the (shock 
frame) Mach number (see equation 13). Calculation of the 
post-shock properties is then simply a matter of evaluat- 
ing the Rankine-Hugoniot jump conditions (equations 8-11). 

4. ) If the mass ratio is different from unity, distribute 
the energy over two shocks. Our simulations suggest that 
the two shocks contribute comparably to the final entropy 
(see Figs. 2, 4, & 8; note also that as the mass ratio 
approaches unity, the double shock model converges to 
the single shock result). To first order, therefore, one can 
assume the shocks each thermalise half of the total energy 
estimated in steps (1) and (2) above (note, however, that 
the results are not very sensitive to exactly how the energy 
is distributed over the two shocks). For the first shock, 
one can calculate the post-shock conditions as in step (3). 
Next, assume the systems adiabatically expand and the 
density drops to approximately 20% of its pre-merger value 
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for the primary system and to some appropriate value for 
the secondary system (see Fig. 23). Using the post-shock 
entropy from the first shock, this new density, and the 
appropriate value for v in (i.e., which corresponds to the 
remaining half of the total thermalisation energy), one can 
calculate the Mach number of the second shock (equation 
13). The final post-shock conditions are then determined as 
usual via the jump conditions. 

5.4 What next? 

These simulations have allowed us to develop a good under- 
standing of how entropy is generated during cluster mergers. 
We have provided an algorithm that encapsulates this phys- 
ical process. The practical application of this is to be able to 
predict the evolution of the entropy profiles of groups and 
clusters as they grow in mass in a CDM universe. In future 
papers, we will consider this in detail. In particular, we will 
focus on the problem of explaining the self-similar growth 
of clusters seen in hydrodynamical simulations. 

The key application of our results is to understand how 
perturbations in the entropy distribution of clusters prop- 
agate through the merging hierarchy. The entropy profile 
is modified both by cooling, which lowers the entropy of 
the system, and non-gravitational heating (e.g., from super- 
novae, AGN outflows, thermal conduction), which raises it. 
In addition to the inclusion of merger rates derived from 
numerical simulations, this aim requires us to validate the 
heating model developed here using simulations of merg- 
ing clusters with "perturbed" initial entropy profiles (i.e., 
where gas does not trace dark matter). Among others, we 
will explore common physically-motivated examples of en- 
tropy modification include shifting and truncating the dis- 
tributions (e.g., Babul et al. 2002; Voit et al. 2002). This 
study is currently underway (McCarthy et al. , in prep.). 

Aside from such practical applications, there is remain- 
ing academic work to be done as well. In the current study 
we have presented a detailed exploration of how entropy is 
generated in merger shock heating events. However, why the 
entropy is generated in this fashion needs further clarifica- 
tion. For example, an interesting result of our study is that 
the gas-only mergers preserve self-similarity. It is well-known 
that dark matter-only simulations also approximately pre- 
serve self-similarity through the hierarchy (e.g., NFW). Why 
it should be that the gas, which exchanges energy with it- 
self through shock heating, and the dark matter, which ex- 
changes energy with itself through phase mixing and violent 
relaxation, both give rise to the same equilibrium state is 
not immediately obvious (see, e.g., Faltenbacher et al. 2006). 
Presumably, this is the result of both the gas and dark mat- 
ter adhering to the virial theorem, but demonstrating this 
explicitly is non-trivial. Another interesting result is that the 
simplest of shock heating models, a single-shock model, does 
not provide an adequate description of mergers characterised 
by large mass ratios. What fundamental factor determines 
how many shocks are necessary to completely thermalise 
the available energy? One possibility is that the system is 
following the course of maximum entropy generation. For 
example, we have found that if the density drop between 
shocks is determined by the amount of kinetic energy re- 
maining to be thermalised (i.e., the leftover kinetic energy 
fixes the degree of adiabatic expansion between shocks), the 



maximum amount of entropy generated almost corresponds 
to the case where the total energy is thermalised equally 
over two shocks. This also appears to be roughly the route 
the simulated systems are following. These and other basic 
matters require further attention. 
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Table 3. Mass Resolution Study. 



Sim. label 






^-gas 


"Mm 








(M e ) 


(Mq) 


lowest res. 


10 4 


1.5 x 10 4 


2.7 x 10 10 


1.7 x 10 11 


low res. 


3.3 x 10 4 


5.1 x 10 4 


8.2 x 10 9 


5.2 x 10 1() 


medium res. 


10 5 


1.5 x 10 5 


2.7 x 10 9 


1.7 x 10 10 


high res. 


3 x 10 s 


4.6 x 10 s 


9.1 x 10 8 


5.8 x 10 9 


highest res. 


10 6 


1.5 x 10 6 


2.7 x 10 8 


1.7 x 10 9 



APPENDIX 

I. MASS RESOLUTION STUDY 

Given that our simulations are non-radiative, the number of 
particles required to accurately capture the evolution of the 
systems should not be particularly stringent. However, to 
be sure we have carried out a mass resolution study for one 
of our simulations. In particular, we have simulated the 1:1 
gas+DM small impact parameter merger at several different 
mass resolutions (see Table 3 for a summary). 

In Figure 26 we show the final entropy profiles of the 
merged system for the five different mass resolution runs. 
Instead of radius, we use integrated gas mass along the ab- 
scissa. Both coordinates have been scaled to the anticipated 
self-similar result (see §3.1). 

Figure 26 shows that resulting profile is remarkably in- 
sensitive to the adopted mass resolution. For example, there 
is only a 20% shift between the resulting entropy profiles of 
the lowest and highest resolution simulations, even though 
the resolution differs by a factor of 100 between the two. 
To strike a balance between speed and accuracy, we adopt 
the characteristics of the medium resolution run for all of 
our other simulations. As indicated by Fig. 26, the medium 
resolution run yields a final entropy profile that differs only 
by a few percent from our highest resolution run. 

We point out that the increased entropy in the low- 
est resolution run is likely due to an underestimate in the 
gas density which, in turn, results in more efficient entropy 
generation in the shocks. However, as one increases the res- 
olution (i.e., particle number), one obtains a more accurate 
density determination and, therefore, a more accurate en- 
tropy jump. In the case of high mass ratio mergers, our 
simulations indicate that most of the energy is thermalised 
in the more massive primary system. Therefore, so long as 
the primary system is well-resolved the results should be ro- 
bust. This likely accounts for the fact that the distribution 
of gas in massive virialised systems formed in non-radiative 
cosmological simulations does not depend much on resolu- 
tion (e.g., Frenk et al. 1999), even though the properties of 
the small systems that merge to form the massive system 
change significantly with resolution. 
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Figure 26. Comparison of final entropy profile for the 1:1 
gas+DM small impact parameter merger simulated with various 
mass resolutions. 
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